Effects of internal noise in mesoscopic chemical systems near Hopf bifurcation

The effects of internal noise in mesoscopic chemical oscillation systems have been studied analytically, in the parameter region close to the deterministic Hopf bifurcation. Starting from chemical Langevin equations, stochastic normal form equations are obtained, governing the evolution of the radius and phase of the stochastic oscillation. By stochastic averaging, the normal form equation can be solved analytically. Stationary distributions of the radius and auto-correlation functions of the phase variable are obtained. It is shown that internal noise can induce oscillation; even no deterministic oscillation exists. The radius of the noise-induced oscillation (NIO) becomes larger when the internal noise increases, but the correlation time becomes shorter. The trade-off between the strength and regularity of the NIO leads to a clear maximum in its signal-to-noise ratio when the internal noise changes, demonstrating the occurrence of internal noise coherent resonance. Since the intensity of the internal noise is inversely proportional to the system size, the phenomenon also indicates the existence of an optimal system size. These theoretical results are applied to a circadian clock system and excellent agreement with the numerical results is obtained.


Introduction
In recent years, the effects of internal molecular noise in small-scale chemical reaction systems, where the number of reactant molecules can be low, have gained much attention [1]- [22]. Though noise may be viewed as nuisance at the first glance, more and more studies show that it can also play counterintuitive constructive roles, especially near some kind of critical points in nonlinear systems. Specifically, for mesoscopic chemical oscillation systems, internal noise can induce stochastic oscillation, in the dynamic parameter region outside but close to the deterministic oscillatory region [3,4]. In addition, the effective signal-to-noise ratio (SNR) of the noise induced oscillation (NIO) shows a clear maximum when the noise intensity changes, which has been well known as internal noise coherent resonance (INCR) [5]- [13]. Since the magnitude of the internal noise is proportional to 1/ √ V , where V denotes the system size, the above phenomenon was also called system size resonance (SSR). Since oscillatory dynamics is of ubiquitous importance in nature, such a behavior may find wide applications in a variety of systems, such as calcium signaling [14], ion-channel gating [15], gene expression [16] and circadian oscillation [17]- [19] taking place in subcellular spaces, as well as catalytic reactions on the surface of nanoparticles or single crystals [20]- [22]. Most of the studies are based on numerical simulations; however, the common mechanism of such behaviors is still unknown to us.
In the present paper, we have developed a theory to understand the NIO and INCR. We notice the fact that they all happen close to the Hopf bifurcation (HB) point of the deterministic system, indicating that some common features of HB must be relevant. Based on the bifurcation theory of vector fields, the behavior of the system near HB can be described by a normal form equation on a two-dimensional (2D) center manifold [23]. Starting from the chemical Langevin equations (CLEs) recently proposed by Gillespie [24], a stochastic normal form equation is obtained near the HB. Thanks to the method of stochastic averaging, the stochastic normal form can be solved analytically, giving the exact expression of the most probable radius r s , auto-correlation time τ c , and SNR of the NIO. When the system size V decreases (the internal noise increases), r s increases monotonically, i.e. the NIO becomes stronger, while τ c decreases monotonically, indicating that the NIO becomes less correlated in time. To a good approximation, the effective SNR equals to (r s τ c ) 2 , which undergoes a clear maximum at an optimal system size V opt . This analysis gives a general picture of the effect of internal noise near 3 the HB: it destroys the regularity of the oscillation induced by itself. The 'trade-off' between the strength and regularity of the NIO leads to the resonance phenomenon.
The paper is organized as follows. Section 2 gives the theory. In section 3, we apply the theory to a circadian clock system. Discussions and conclusions are drawn in section 4.

The CLE
Consider homogeneous chemical reactions involving N species and M reaction channels taking place in a small volume V , Each reaction channel ρ changes the integer numbers of X i ∈ N of molecules of the species i by an amount equal to the stochiometric coefficient The transition rate associated with each reaction reads, It is now generally accepted that internal fluctuations, resulting from stochasticity of the discrete reaction events, become considerable in small systems. Generally, one may view the reaction as a Markovian stochastic process, and the system's dynamics can be described by a master equation ruling the evolution of the probability P({X i }, t) for the system in state {X i } at time t: Very recently, Gillespie argued that if some kind of 'macro-infinitesimal time scale' exists in the system, the master equation can also be approximated by a CLE [24]. Generally, if the system size is not too small, such a condition is expected to be fulfilled and the CLE works, at least in a qualitative manner. In practice, the CLE is much faster for computer simulation, and based on some previous studies [3,4], the CLE showed rather good agreement with exact simulation methods of the master equation. For the reaction network (1), the CLE reads, Here ξ ρ = 1,2,...,M (t) stand for independent Gaussian white noise with ξ ρ (t) = 0 and ξ ρ (t)ξ ρ (t ) = δ ρρ δ(t − t ). If the molecule population is not too small, one may approximately write the transition rate as 4 where x i = X i /V is the concentration of the ith species. Hence the CLE (2) changes to, It is clear that in the thermodynamic limit V → ∞, one recovers the deterministic equation for reactions (1), Therefore, the second term on the right-hand side of equation (7) stands for the 'internal noise', the magnitude of which is proportional to 1/ √ V . Also note that the internal noise is coupled into the system's dynamics in a multiplicative manner, and each reaction channel involves an independent Gaussian white noise ξ ρ (t). Equation (7) will be the starting point of our following analysis.
One should note that there are certain limitations on the use of the Langevin approach. As stated above, Gillespie emphasized the existence of a 'macro-infinitesimal' time scale for the validity of CLE. This condition cannot always be true in real systems such that the predictions of the CLE could show apparent discrepancies with the experiments or simulation results. For instance, in the case of the propagation of chemical fronts into an unstable state, the Langevin method predicts that the correction to the front speed is positive and scales as V −1/3 [25], while microscopic simulation and master equation approaches demonstrated that the correction is negative and scales as (lnV ) −2 [26]. One should also note that Gillespie's version of the Langevin equation for mesoscopic chemical reactions was not the first one, as also stated in Gillespie's paper [24], and it is associated with the second order truncation of the Kramers-Moyal expansion for the reaction scheme [27]. Despite these limitations, however, the agreements with numerical simulations in the present work indicate that the CLE is valid to study the phenomenon of NIO and INCR, even in some quantitative manner.

Stochastic normal form
As stated above, we are interested in the effects of noise for parameter regions close to the supercritical HB of the deterministic system (8). We denote the control parameter by µ and the HB by µ c . Without losing generality, we assume that oscillation happens for µ > µ c . According to the Hopf theorem [28], the Jacobi matrix (J) i j = ∂ F i ({x i })/∂ x j has a pair of conjugate eigenvalues λ ± = α ± iω for µ µ c , with α < 0(>0) for µ < µ c (> µ c ). The other N − 2 eigenvalues of J all have strictly negative real parts with absolute values considerably larger than 0. Hence near the HB, all those N − 2 modes will relax much faster, and the system's dynamics will be dominated by the slow motion on a 2D manifold tangent to the subspace spanned by the eigenvectors of λ ± . Then by variable transformation, one can define a complex amplitude Z , which obeys the following 'normal form' equation, where C r and C i are constants determined by the nonlinear terms in F i = 1,...,N ({x i }). We consider that the HB are supercritical, such that the oscillation is stable and C r < 0.
We now turn to the CLE (7) which includes the internal noise terms. Due to the presence of noise, one should be careful during the variable transformation. One notes that to be consistent with the master equation, equation (7) must be interpreted in the Ito manner. To facilitate the variable transformation, one should first transform it into a Stratonovich form, which allows for normal calculus [27]: where • stands for the Stratonovich interpretation, and It is clear that the difference resulting from this transformation is of the order of 1/V . Although the physical concept should be clarified here, qualitatively, this small correction can be neglected. We thus will simply set F = F in equation (10) in the following analysis.
For the deterministic part of equation (10), we can then follow the standard procedure to get the normal form. During this process, the noise terms also undergo a linear transformation. After some tedious but not difficult manipulation, one finally obtains the following 'stochastic normal form' [29] for Z , Herein, coefficientsν 1ρ andν 2ρ are the elements of the transformed stochiometric matrix (ν) iρ = (T −1 ν) iρ . Writing Z = r e iθ , where r and θ can be viewed as the amplitude and phase of the oscillation respectively, we have, where χ rρ = (ν 1ρ cos θ +ν 2ρ sin θ ) √ w ρ and χ θρ = (−ν 1ρ sin θ +ν 2ρ cos θ) √ w ρ /r . In the vicinity of the supercritical HB (|α| 1), the time scale for r and θ can be separated. This fact makes it possible to use the 'stochastic averaging' procedure [30], which can approximate the system as Markov processes in the long time limit. Consequently, one can approximate equation (13) by the following Ito stochastic differential equations: resulting from the coupling between r and θ . ε 2 r = ρ 2π 0 χ 2 rρ dθ/2π and ε 2 θ /r 2 = ρ 2π 0 χ 2 θρ dθ/2π are the averaged noise intensities associated with the two 'new' independent Gaussian white noises ξ r and ξ θ . Note that we can expand the reaction rates w ρ in the following way (because the state variable x is related to r and θ via the linear transformation), If the reactions are 'elementary' processes and assuming the validity of the mass-action law, n will be the maximal 'order' of the reactions which is usually not greater than 3. Sometimes the reactions in equation (1) are combined from many elementary steps via quasi-steady state approximation and the mass-action law does not hold, n could be arbitrary. However, since r is small near the HB, only a few leading terms contribute to w ρ . Substituting equation (16) into (15), we can easily find that K (θ) is zero and only those coefficients with even k + l have nonzero contributions to K (r ), ε 2 r and ε 2 θ . The latter ones all have the form ε 2 + (kl) γ (kl) r k+l , where the summation runs over terms with even k + l only, γ (kl) are constants determined by w (kl) ρ and For small internal noise level, r 2 1 and it is a good approximation to neglect the terms with k + l 2. Therefore, equations (14) finally reduce to a rather simplified form, We may call this equation the 'stochastic averaged normal form equation (SANFE)' of the CLE (7).

NIO and coherent resonance
The SANFE already shows new features of the system's dynamics. In the first equation, a 'new' deterministic term ε 2 /2V r appears, which vanishes in the macroscopic limit V → ∞. Even for α < 0, where no oscillation happens in the macroscopic system, one can still have nonzero solution for αr + C r r 3 + ε 2 /2V r = 0, which is Therefore, some kind of oscillation happens that is 'induced' by the internal noise, say, NIO. In addition to a finite system size V , it is shown that a nonzero w (00) ρ is also necessary for the occurrence of NIO. This is always true for chemical reactions with nonzero steady state concentrations x s . Therefore, NIO is a 'universal' behavior for mesoscopic chemical systems staying outside, but close to, the deterministic supercritical HB.
From the SANFE, we can readily write down the Fokker-Planck equation for the probability distribution function of r , By setting ∂ t P s (r, t) = 0, we can get the stationary distribution function where C 0 is a normalization constant. p s (r ) has a maximum at r = r s , hence in the stationary state, the system will most probably stay around the noise-induced limit cycle. Approximately, we may then replace r 2 and r −1 in the second equation of (18) by r 2 s and r −1 s respectively. θ is then Gaussian distributed with mean and variance, Using e iθ = e i θ −[ θ(t) 2 − θ(t) 2 ]/2 for Gaussian random variable [27], and after some straightforward mathematical manipulations, we can get the auto-correlation function (ACF) for cos θ as following, The ACF is a typical damped oscillation, with frequency given by ω 1 and correlation time given by τ c . As for the ACF C(τ ) of the state variable y 1 r cos θ, we can simply multiply C θ (τ ) by r 2 s , to keep the leading term (one can also keep the terms of higher order, but that only lead to minor quantitative correction and is not necessary). By Fourier transformation of C(τ ), we can get the power spectrum density (PSD), which has a Lorenzian-like form. Clearly, the PSD has a peak at ω = ω 1 , whose height H and half-height width ω are given by In the literature, the performance of the NIO is often characterized by the effective SNR, which is defined as the peak height divided by the half width [4], hence, It is easy to verify that ∂r s /∂ V < 0 and ∂τ c /∂ V > 0. Therefore, when the internal noise level increases (system size V goes small), the amplitude of the NIO becomes larger, while the correlation time decreases. In other words, the strength of the NIO is enhanced, but the regularity is reduced. One may expect that for some optimal system size V , the 'trade-off' between the 8 strength and regularity of the NIO gives a maximal SNR. By ∂(SNR)/∂ V = 0, the optimal system size reads, It is interesting to note that at this optimal size, the radius r s of the NIO fulfill the equality (αr s = C r r 3 s ) V =V opt i.e. the linear and nonlinear terms contribute equally to the evolution of r in the original deterministic normal form equation.
We conclude that the above procedures apply to any reaction system with supercritical HB. One can just follow some easy steps to use the theory. To begin, write down the deterministic dynamics equation and find the HB by use of any convenient bifurcation software or just computer simulation, put your control parameter close to the HB but in the steady state region. Then, calculate the eigenvectors of the Jacobi matrix, which can be easily done by Matlab or self-written codes. After that, construct the transformation matrix T = (Re u + , −Im u + , u 3 , . . . , u N ), and calculate the matrixν = T −1 ν; this gives the valuesν 1ρ and ν 2ρ for ρ = 1, 2, . . . , M. Finally, equations (21)- (28) give the final results.
Before ending this section, we would like to emphasize that the stochastic normal form and averaging procedure have been used to study many stochastic dynamic systems in the literature ( [29] and references therein). Nevertheless, no studies so far have considered the behavior of NIO and INCR. An important point is that our study is based on a real chemical reaction system, rather than a general mathematical model. Recalling the procedures above, we know that the term ε 2 /2V r plays the key role for NIO and INCR. Based on the expression ε 2 = ρ (ν 2 1ρ +ν 2 2ρ )w (00) ρ /2, this term will be exactly zero if the equilibrium points x s are zero. For chemical reactions, x s stand for steady state concentrations of some species and cannot be zero. This is why NIO and CR were not reported in some literature regarding stochastic perturbed HB like [31], because a zero fixed point was assumed 'without loss of generality' therein and the only consequence of noise was a slight shift of the Hopf point. Therefore, our analysis reveals that NIO and INCR are features of real chemical oscillating systems rather than of an arbitrary mathematical model.

Application to circadian clock system
In this section, we will apply our theory to a circadian clock system. Most living organisms use circadian clocks to keep an internal sense of daily time and adapt their behavior accordingly. It is now known that the circadian clock system is regulated by a gene network on the molecular level, such that internal noise must be considered. Actually, many studies have been made of this issue, but most works so far have assumed that the internal noise is destructive and they mainly focus on the robustness or resistance of circadian oscillations to such internal noise. Back in 1963, Goodwin had considered the effect of intrinsic noise in genetic oscillator models and he noted that a minimum number of molecules was required for the clock robustness [32]. In [33], Gaspard performed a theoretical study on the robustness of mesoscopic chemical clocks, and he found that a minimum number of molecules is required for the mesoscopic oscillations to remain correlated in time. Barkai and Leibler argued that the sensitivity to internal noise and the robustness to such uncertainties were probably decisive factors in the evolution of circadian clocks, and should be reflected in the underlying oscillation mechanism [34]- [36]. In a recent paper [3], on the other hand, we have focused on the 'constructive role' of internal noise and Translation of P C out of the nucleus demonstrated that INCR phenomena, as predicted theoretically in the previous section, did happen in this system. Therefore, the main purpose here is to compare the simulation results with the theory using this particular system. The minimal model used in the present paper incorporates the transcription of the gene (G) involved in the biochemical clock and transport of the mRNA (R) into the cytosol where it is translated into clock proteins (P C ) and degraded. The protein can also transport into the nucleus (P N ) where it exerts a negative regulation on the expression of its gene. Such negative regulation forms the core mechanism of the oscillation behavior [18]. Based on this mechanism, there are six reaction channels, as listed in table 1.
Following the procedures in the last section, we first need to locate the HB for the deterministic system. By numerical simulation, we can easily locate the approximate region where the bifurcation happens. Then, we carefully change the control parameter in this region, calculate the fixed point x s = (x 1s , x 2s , x 3s ) and the eigenvalues λ ± = α ± iω of the Jacobian matrix. The HB can then be exactly located when α passes zero. For the parameters chosen, we find the HB at ν HB s = 0.25725 ± 0.00002, above which oscillation happens. NIO and coherent resonance are expected to happen to the left-hand side of ν HB s . For a given v s , the fixed point x s , eigenvalues (λ ± = α ± iω, λ 3 ) and the corresponding eigenvectors (u ± , u 3 ) of the Jacobian matrix are calculated. Via normalization, we can set the first nonzero component of u + to 1. Write the eigenvectors as u ± = (1a ± ibc ± id) , u 3 = (e f g) (here the superscript denotes vector transpose), the transformation matrix reads One can then calculate the matrix T −1 , and byν = T −1 ν the coefficientsν 1ρ andν 2ρ for (ρ = 1, 2, . . . , 6) are obtained. C r and C i can be calculated numerically, and the result is C r −0.3474 and C i 0.5772. From these numerical values, the stationary probability distribution P s (r ), ACF C θ (τ ), the correlation time τ c , and SNR can be obtained theoretically from equations (21)- (23) and (27). These will be compared with direct numerical simulations of the CLE (7).
Here are some details of the numerical simulations. For the CLE (7), the Ito interpretation is used and numerical integrations are performed according to standard methods of stochastic calculus [27]. We choose t = 0.001 in our simulation. After a long transient time, r (t) is calculated as r (t) = y 2 1 (t) + y 2 2 (t), where the vector y = (y 1 , y 2 , y 3 ) are transformed from x = (x 1 , x 2 , x 3 ) via y = T −1 (x − x s ). Accordingly, cos θ(t) = y 1 (t)/r (t). The probability distribution of r (t) is calculated over a long enough time period. PSD are calculated from the time series of y 1 (t) with 16 384 data points filtered by a Welch window. We perform a nearestaveraging smoothing over 50 points on the PSD curves to estimate the SNR values.
In figure 1(a), the comparison of the stationary distribution of r (t) between the theory and numerical simulation, for ν s = 0.257 (α = −0.00038), is shown. Reasonable agreements are observed, especially for small noise levels (large system size). With the increment of internal noise level, the distribution becomes wider, and the most probable value of the NIO r s becomes larger. The dependence of r s on system size V is depicted in figure 1(b), where simulation and theory show excellent agreement. In figure 2(a), the comparison is shown between the ACF of cos θ (t) obtained from simulation and equation (23), for v s = 0.250 (α = −0.0126), and they also show reasonable agreement. By peak fitting the ACF from the numerical simulation, we can estimate the correlation time as a function of the system size, as shown in figure 2(b). Also shown in figure 2(b) is H = r 2 s τ c , representing peak height in the PSD indicated by equation (25). Good agreement between the theoretical predictions and simulation data can be observed. Finally, the dependence of SNR on the system size is shown in figure 3(a), by numerical calculation, and figure 3(b), by equation (27), for ν s = 0.257 (α = −0.00038), ν s = 0.255(α = −0.00355) and ν s = 0.250 (α = −0.0126), respectively. We note that the qualitative features agree with each other, i.e. INCR appears and the optimal size V opt and the maximal SNR both become larger when the distance from the HB decreases. Quantitatively, there are   (27)). still discrepancies regarding the exact location of the optimal system size V opt and the maximal value of SNR (note that the absolute values of the SNR do not make sense because the SNR are in arbitrary units). The main reason is that correct numerical estimation of the PSD and henceforth SNR of the noisy data is difficult. As a whole, we can draw the conclusion that our theory reproduced well all the numerical results.
Before ending this section, it is worth emphasizing that although being widely used in the literature, the so-called SNR defined as 'peak height divided by the width' is somewhat misleading. Actually, it is not a ratio of signal to noise, but rather a measure of the performance of the NIO. As already pointed out in the last section, the maximum in the SNR simply represents a type of 'balance' between the regularity and strength of the NIO. According to 12 the analysis there, both the strength (denoted by the peak height H in the power spectrum or the oscillation amplitude r s ) and regularity (denoted by the correlation time τ c or equivalently 'width of the peak ∆ω divided by the mean frequency ω 1 ) of the NIO change monotonically with the system size V without any 'resonance'. For too small internal noise (large V ), r s is very small and NIO can hardly be observed although the calculated τ c could be large. On the other hand, if internal noise is too large (small V ), NIO can hardly be distinguished from 'random noise' because τ c is too small. Only for intermediate noise levels (optimal V ), may the observed NIO have both considerable strength and long correlation time, corresponding to INCR. In figure 1(b) and 2(b), we have shown that r s , H , and τ c all change monotonically with V , in contrast to the SNR drawn in figure 3, where clear maxima exist.

Discussions and conclusions
In the literature, the effects of noise on mesoscopic chemical oscillations have also been studied analytically by others. Vance and Ross [37] have investigated the fluctuations near the limit cycle using a master equation approach. They reduced the master equation to a Hamiltonian-Jacobi equation in the large size limit, and constructed approximate time-dependent as well as stationary solutions to the master equation in the vicinity to the limit cycle. Gaspard et al have studied the correlation time of mesoscopic chemical clocks also using the master equation [33,38], and an estimation was obtained for the minimum number of molecules required for the chemical oscillations to remain correlated in time. However, their analysis mainly focused on the parameter region where deterministic oscillations already exist, and internal noise was considered to be destructive by inducing phase diffusion. It is worth noting that their analysis may fail when bifurcation happens, as pointed by Gaspard, while our theory only applies to the vicinity of the bifurcation point. Therefore, combining those studies and the present work, one may be able to have a deeper understanding of the effects of internal noise in mesoscopic chemical oscillating systems.
To conclude, we have developed a theory for NIO and INCR in mesoscopic chemical oscillation systems, by using a stochastic normal form equation and stochastic averaging procedure. The analysis is quite general, and can be applied to any chemical reaction system with supercritical HB. We apply the theory to a circadian clock system, and the theory and numerical simulation show reasonable agreement. Since chemical oscillations are of ubiquitous importance in nature, especially in living systems, the analysis in the present paper could find applications in many other systems.