# On the dephasing of genetic oscillators

See allHide authors and affiliations

Contributed by Peter G. Wolynes, December 18, 2013 (sent for review November 20, 2013)

## Significance

Genetic oscillators are ubiquitous regulatory motifs in the molecular control circuits of living cells. Prominent examples include the cell cycle and cellular signaling. There are two primary sources of noise in these oscillators: the binary character of gene states and the low copy numbers of proteins. This molecular noise induces dephasing in oscillators in analogy to the way imperfections in clocks lead to disparities in time. We can study how these two distinct sources contribute to dephasing and how well it can be approximated by adding noise to a macroscopic model. Using the *NFκB*/*IκB* oscillator as an example, we find that gene noise leads to significant deviations from the often-used phenomenological models.

## Abstract

The digital nature of genes combined with the associated low copy numbers of proteins regulating them is a significant source of stochasticity, which affects the phase of biochemical oscillations. We show that unlike ordinary chemical oscillators, the dichotomic molecular noise of gene state switching in gene oscillators affects the stochastic dephasing in a way that may not always be captured by phenomenological limit cycle-based models. Through simulations of a realistic model of the *NFκB*/*IκB* network, we also illustrate the dephasing phenomena that are important for reconciling single-cell and population-based experiments on gene oscillators.

Cyclic rhythms are a common feature of many self-organized systems (1) manifesting themselves in myriad forms in biology, ranging from subcellular biochemical oscillations to cell division and on to the familiar predator–prey cycles of ecology. An oscillatory response of a gene regulatory circuit, whether transient or self-sustained, can provide several advantages over a temporally monotonic response (2, 3). The ability of copies of a system to synchronize may lead to dramatic noise reduction and greater precision for timing in assemblies. On the subcellular level, rhythmic dynamics span time scales from a few seconds, as in the calcium oscillations, to days, as in the circadian rhythms, or years for cicada cycles (1, 4⇓–6). Ultradian genetic oscillations, which take place on an intermediate scale from minutes to a few hours, are medically important. A singularly important case is the Nuclear Factor Kappa B (*NFκB*) gene network, which organizes the mammalian cell’s response to various types of external stress and plays a role in regulating inflammation levels in populations of cells. The response of the *NFκB* circuit to continuous external stimulation has been studied by Hoffmann et al. (7) who observed damped oscillatory dynamics for *NFκB*, which has been linked to the presence or absence of particular isoforms of the inhibitor *IκB*. On the other hand, experiments carried out on individual cells have detected more sustained *NFκB*/*IκB* oscillations that either are completely self-sustained (8, 9) or damp at a much slower rate (9) than found for the population, depending on the duration of external stimulation. It follows that some type of averaging takes place, but the physical mechanisms and the stochastic aspects of this population averaging are not fully understood. Many aspects of the *NFκB* oscillatory dynamics may be rationalized using deterministic mass action rate equations (10, 11), but how stochastic self-sustained oscillations average out at a cell population level remains unclear.

In this work, we provide a conceptual framework for understanding stochastic averaging as a result of “dephasing” of genetic oscillators. By dephasing, one essentially means the loss of common phase or coherence of oscillations in populations of *mRNA* or protein by-products of gene activation, caused by the stochastic molecular events in the course of an oscillator’s operation that occur at different times in different cells. When the phases of genetic oscillators are not controlled by external coupling to other cells or by exogenous signals, the random character of molecular events induces phase diffusion. In the absence of synchronizing forces, an ensemble of genetic oscillators in the long time limit inevitably will be dephased completely as part of the entropic drive to randomize the phase distribution across the ensemble. In an extreme deterministic limit of ”newtonian” genetic oscillators, there would be no dephasing and the oscillators would preserve the memories of their phases forever. The oscillators of the cell, on the other hand, are driven by inherently fluctuating molecular entities and will get out of phase and eventually lose memory of their initial phase in a rather modest time (12, 13).

We explore a particularly simple yet realistic model of the *NFκB*/*IκB* circuit (Fig. 1) and demonstrate how self-sustained stochastic single-cell oscillatory dynamics yield damped oscillations at the population level, as observed in the experiments of Hoffmann et al. (7). Another related but more fundamental question is how the single-molecule nature of the gene contributes to the stochasticity of gene oscillator networks. In our model, we explicitly account for the highly non-Gaussian noise that comes from a single gene switching on and off and investigate how the time scale of gene-state fluctuation affects the noisiness of the oscillations.

## Discrete-State Markov Models of Dephasing

Genetic oscillators are chemically driven mesoscale systems operating with a finite number of proteins and a handful of gene states. A satisfactory description for the dynamics of gene oscillators is given by a master equation, which captures the underlying discreteness and the probabilistic nature of molecular transactions. In such a description, the state of the system *S* is given by specifying the number of each type of protein together with the occupation of states of the genes. Once transition rates *W*_{ij} between the pairs of states *i*, *j* ∈ *S* are assigned, the probability ket can be described starting from a given initial condition . Assuming Markovian dynamics for the transitions, we haveThe elements of the stochastic rate matrix **W**, are the transition rates from states *j* to *i*— for *i* ≠ *j*—and the net escape rates from states *i—*. The probability of a microscopic state, *z*, then is given by . Owing to time-translation invariance, we can write the solution using the eigenvectors and eigenvalues of the rate matrix **W**. Assuming the eigenvalues are not degenerate, the general solution is . The conditional probabilities for states *z* are then given by . The Perron–Frobenius theorem requires the existence of at least one purely real eigenvalue with real part zero, whereas the rest of the eigenvalues must have strictly negative real parts . For systems with broken detailed balance, the rate matrix is non-Hermitian (14), which may produce complex eigenvalues corresponding to oscillatory motion along *N*_{c} states for which ratios of left and right fluxes do not balance (14, 15). Circular fluxes do not average out when system size is scaled, which leads to limit cycles on macroscopic scales (16). A convenient way to quantify temporal relaxation of an ensemble of oscillators to the steady state is via autocorrelation functions of the form . The steady state probability is obtained by simply taking the limit of very long time: . The last equality follows from the fact that the left eigenvector corresponding to the stationary state is unity, . Plugging in the expressions for the probabilities in the expression for and integrating over variables *z*_{τ} and *z*_{0} we obtainwhere the α’s and β’s are defined as , and . Because **W** has real elements, the eigenvalues and eigenvectors come in complex conjugate pairs, i.e., and . Pairs with smaller correspond to the slow evolving modes, which for the oscillating system, also should have a large enough frequency component so that the system goes through at least one cycle before decaying to the steady state. In other words, one expects the mode(s) with the highest and lowest values to dominate oscillatory dynamics. If we assume there is a significant spectral gap between the first eigenvalue and the remaining ones (∀ *i*, and ), we can subtract the constant stationary fluctuation term, , to obtain a fairly simple expression for the autocorrelation function:where the and are constants set by the initial condition. The slowest modes, according to our assumption, are linked to the motion along the cycle, whereas the faster modes come from fluctuations longitudinal to the cycle. The apparent similarity of the form of correlation function **3** to that of more phenomenological approaches (see the next section, Eq. **11**) provides a justification for using such models, which are based on rather restrictive approximations. In the master equation formalism, one needs only invoke the much milder assumption that there is a sufficient spectral gap in the eigenvalues. The last assumption can be tested readily on minimalist models of genetic oscillators. To illustrate the idea, we devised a simple model of a *NFκB* gene oscillator small enough to allow an exact numerical solution of the master equation but complex enough to include all the essential features of gene oscillators. The dynamical variables describing the circuit are the gene state *s*, *NFκB*, number *x*, *IκB*, number *y*, and *mRNA* number *z*. Because there are no reactions that produce or consume *NFκB*, its total amount is conserved (*n* = *const*), which allows full definition of the system through only four variables. A schematic representation of the circuit is depicted in Fig. 1. The rate coefficients of reactions are chosen to ensure that at the macroscopic level, the corresponding mean-field ordinary differential equation (ODE)-based model exhibits limit cycle behavior (see Table S1). Stochastic simulations of a more complete model with variables accounting for cytoplasmic and nuclear concentrations of *NFκB*/*IκB* and with experimentally determined rate coefficients are presented in the last section. The master equation of our simple model has the following explicit form:The *k*’s are rate coefficients and Ω is the dimensionless volume or the size of the system. The terms inside the brackets correspond to the diagonal, and the terms outside the brackets account for various off-diagonal elements of the rate matrix **W**. Numerically, the problem boils down to obtaining a spectral decomposition of a square non-Hermitian matrix with *s* × *x* × *y* × *z* rows (columns). Whereas the variables *s* and *x* have a priori clear upper limits ( and ), no such limit exists for *y* and *z* ( and ). Therefore, to obtain meaningful solutions, one must place boundary conditions and choose the upper values carefully so that transition probabilities beyond the boundaries will have a marginal impact on the qualitative features of the spectrum. We have carried out numerical computations on grids with Ω = 1, *n* = 2 − 10, and , obtaining consistent solutions in all cases. Physically, such an approximation may be justified by realizing that having fewer copies of *NFκB* makes states with larger numbers of *IκB* and *mRNA* molecules less probable; hence, probabilities fall off sufficiently rapidly near the imposed upper boundaries.

The computed spectrum (Fig. 2) reveals the existence of a distinct rotational mode separated from the rest by a relatively large spectral gap. This mode also corresponds to the real eigenvalue smallest in magnitude (or least negative) and therefore the slowest component of the dynamics. All other rotational modes have higher decay rates and, as a result, contribute very little to the observed oscillations (Fig. S1). The spectral gap would be expected to grow linearly proportional to the size of the system Ω based on the WKB expansion-like theories (17, 18) and stochastic simulations (18, 19) of other oscillating biochemical networks. In the following sections, we show that in the case of gene networks, besides macroscopic size, the spectral gap also is a function of the time scale for gene-state fluctuations. For such gene networks, the precise relationship between the size of the system and the spectral gap has yet to be worked out; however, qualitatively one can demonstrate the trends using simplified toy models. As a prototype for such cyclic processes, we consider a cyclic graph composed of discrete states connected by one-way transitions, e.g., . For the cases of equal rates , the rate matrix of such a directed cycle admits analytical solution, yielding the eigenvalues, . Going to the limit of large *N*, we have and . Therefore, as expected, the spectral gap in this model grows linearly proportional to the number of states in the system . The exact solution for the single-gene master equation has been worked out for a few special cases (20). These studies show that slow operator-state fluctuations broaden and skew the product number distributions. To see how slow gene fluctuations accelerate the relaxation to steady state, one must have a model with at least two time scales, one accounting for cyclic flux and another for the gene-state fluctuations coupled to it. As such, we consider a minimalistic toy model of a directed cycle with a unit attached to an edge (Fig. 3).

The role of the edge state is to mimic the binary nature of gene states. This dichotomic variable moderates the cyclic dynamics to an extent depending on the switching time scale. By increasing the rate *G*, the real part of the smallest nontrivial eigenvalue increases, eventually saturating for high switching rates, implying that dephasing dynamics become *G* independent. The spectral solution of Eq. **4** shows the same trend for the values of *k*_{on} and *k*_{off} corresponding to the limit cycle solution in the macroscopic limit. The small system sizes for which we can afford to solve the master equation (Ω = 1) severely limit such explorations, because the oscillators are quite noisy to begin with and dephasing takes place during very short times (see Fig. S1). Fortunately, the study of dephasing in the large () and biologically more relevant limits may be carried out via kinetic Monte Carlo simulations, which are elaborated in the last section.

## Phenomenological Model of Dephasing

In the near-deterministic limit, the complete master equation treatment of gene networks becomes computationally intractable because of the steep rise of the dimensionality of the rate matrix. Fortunately, self-sustained oscillators in the deterministic limit may be described with an autonomous system of equations with limit cycle attractors (21). This opens an avenue for including the noise in a top-down fashion (22, 23). For such attractors, trajectories rapidly relax toward a limit cycle, but once on that orbit, the phase undergoes “free diffusion” driven by the underlying fluctuations. Let us first discuss the deterministic model of self-sustained oscillations that are born via a supercritical Andronov–Hopf (AH) bifurcation (21). Classical examples (1, 4, 5) of chemical oscillations created via AH bifurcation include the Belousov–Zhabotinsky reaction, the Brusselator model, the Selkov model of glycolysis, and many others. The normal form of the AH bifurcation may be written (21)where *z* is the representation of the limit cycle in the complex plane in the vicinity of the bifurcation. In the context of gene circuits, *z* is the dynamical variable that describes the oscillatory cycle formed by any two components in the feedback loop. Such a simple description is a consequence of planarity of the attractor in the phase space. The parameters *p*_{1} and *p*_{3} are complex expansion coefficients , Λ is the bifurcation control parameter, and Λ_{c} is its critical value. Transforming the variables to polar form, one obtains the normal forms for amplitude , and phase . Because of the attractor nature of the limit cycle, the amplitude equilibrates on a faster time scale, whereas the phase evolves freely. Therefore, we may put *r* at its equilibrium value so that the long time behavior of phase will beDeterministic limit cycle dynamics may be described as a combination of an oscillation with a constant angular frequency , and amplitude *r* = *r*_{0}. The accuracy of the deterministic picture deteriorates as we scale down the size of our system so that at some point, we must account for the fluctuating molecular nature of the oscillator (24, 25). As a first-order approximation, we may describe the dynamics of phase and amplitude as being driven by a Gaussian white noise ,

where is the noise-induced shift of the phase and *D*_{r} and *D*_{ϕ} are the diffusion coefficients for amplitude and phase. The λ and α are parameters that depend on the rate coefficients in the gene network. For the limit cycle to be stable with stationary amplitude , the parameters must be strictly positive α, λ > 0. It is convenient to work with the corresponding Fokker–Plank equation . Writing , the Fokker–Planck operator readsBecause the dynamics of the amplitude are decoupled from phase diffusion, the *P*(*r*, *ϕ*, *t*) becomes simply a product of two individual distributions . Because we also assume that amplitude dynamics take place on a time scale faster than the phase fluctuates, we may use the steady-state probability distribution satisfying the equation . The phase distribution is seen to be that of a simple Wiener process and is given by the standard form of Green’s function for the diffusion equation with modulo 2π. Taking the time evolution from a specified initial condition of phase , one obtains the following form for the probability distribution:where the *N*^{−1}(*t*) is a normalization factor. At this point, it becomes straightforward to compute the correlation function by averaging the time-lagged product of observables, . As phase diffusion is decoupled from amplitude fluctuations, we can perform the averages separately and take the product at the end. The radial part yields . Numerical evaluation of the radial part, assuming some reasonable values for the parameters, shows that , implying that noise increases the effective length of the oscillation amplitude. For the phase, it is more instructive to derive the correlation function by first writing the phase, . In the complex exponential representation, after taking the average over realizations of noise using the Gaussian white noise assumption, one finds . Combining the last two expressions for amplitude and phase, we obtain the final expression for the correlation function:The correlation function is a damped cosine oscillating at an average stochastic frequency ω with a constant average amplitude. In this mesoscopic description of dephasing, the damping time scale is set by the noise intensity *D*_{ϕ} of the phase variable. It represents a “virtual” damping, because it results from dephasing of the trajectories due to stochastic fluctuations, whereas individual trajectories themselves would appear to continue to oscillate. Finally, to connect the phenomenological picture with the master equation-based description of the previous section, we express the solution of the phase diffusion equation via eigensolutions of the equation , where and the eigenvalues and eigenfunctions are and , respectively. From the expression of eigenvalues, we find the spectral gap to be , where the dephasing coefficient is inversely proportional to the size of the system, (17), making the oscillations completely dominated by the main rotational mode. The expression for the correlation function takes the following form:where *C*_{n} are constant parameters set by the initial conditions. In the near-macroscopic limit , one recovers Eq. **11**.

## Stochastic Simulations of Genetic Oscillator Models

At a simplified level, the core of the NFκB gene network involves a negative feedback (Fig. 1) by the inhibitors *IκB*, which suppress the binding of the NFκB to the gene that activates the synthesis of the IκB. In eukaryotic cells, the activation of *NFκB* is initiated by external cellular stimulation, which triggers a chain of reactions leading to the synthesis of IKK proteins. The IKK phosphorylates the complex of *NFκB* with *IκB*, marking the *IκB* for degradation and thereby activating the *NFκB*. Once freed from *IκB*, the *NFκB* diffuses to the cell nucleus and activates hundreds of genes, among which is also the *IκB* gene, which down-regulates the activity of *NFκB* to ensure the attainment of a homeostatic state once the external stimulus is over. The cyclic suppression of *IκB* synthesis by *IκB* itself creates a closed time-delayed loop, resulting in oscillatory behavior. Our model of the network is based on the scheme proposed by Sneppen and coworkers (10, 11), who explored the oscillatory behavior of *NFκB*/*IκB* using a system of deterministic ODEs. A key difference between their model and the present one is that we include the digital nature of gene activation (*OFF* + *NFκB* → *ON*; Fig. 1) as opposed to modeling transcription using pseudo–first-order mass-action kinetics with Hill coefficients. Although sharing the main qualitative features, the model we use for stochastic simulations is more sophisticated than the toy model of an *NFκB* oscillator used to solve the master equation. Most importantly, this more sophisticated model accounts for the cytoplasmic and nuclear concentration of all components and uses rate coefficients extracted from the experiments (see *SI Text* and Table S2 for details). The stochastic simulations are set up to mimic the experiment of Hoffmann et al. (7) in which the cells are under constant exposure to external stimuli, which results in steady production of *IKK* and keeps the single-cell oscillations undamped, whereas the cell population average decays over a finite period (9).

We stochastically simulate the system of reactions via a kinetic Monte Carlo algorithm (26) (see *SI Text* for details). Although both deterministic (see Fig. S2) and stochastic levels of treatment show sustained oscillations in the populations of *NFκB*/*IκB*, stochastic simulations have gene noise-induced dephasing, which damps the oscillations at an ensemble level (Fig. 4 and Fig. S3). Averages are taken over the independent stochastic trajectories initiated from the same state. The latter is done to dissect the effects of gene noise on the dephasing by eliminating the extrinsic noise. In contrast to the master equation studies, here we use a much larger system size (Ω ∼ 10^{3}), which is a better approximation of the biochemical environment in which *NFκB IκB* gene oscillators operate in mammalian cells. Gene noise is amplified by reducing binding/unbinding rates, which leads to faster dephasing in the ensemble of oscillators. In fact, theoretically one may achieve the same levels of dephasing observed in the experiments of Hoffmann et al. (7) by simply decreasing the binding (*k*_{on}) and unbinding (*k*_{off}) rates while maintaining the ratio equal to its known (27) in vitro value , as is done in our simulations. This provides an alternative explanation for the damping of the oscillations from that suggested by the deterministic models, which would indicate that the experimentally observed damping must arise through the action of the other members of the *IκB* family via independent reaction pathways.

We see that the single-molecule nature of the gene in our model turns out to be quite crucial for describing the dephasing of genetic oscillators, because the time scale of operator binding and release has a significant impact on the dephasing times τ_{ϕ}. In the language of the discrete Markov and phenomenological models, slowing of operator binding/release transitions boosts the dichotomic gene noise, which makes the real part of the eigenvalue corresponding to the oscillatory mode more negative and increases the apparent phase diffusion (Fig. 4 and Figs. S1, S3, and S4), leading to faster damping of the oscillations.

To better appreciate the nature of the single-gene dynamics in the stochastic dephasing and to test the phenomenological limit cycle-based models, we compute 2D stationary probability distributions by varying the system size Ω (Fig. 5) for the cases of fast (adiabatic) and slow (nonadiabatic) gene binding/unbinding (28) and also in the absence of dichotomic gene noise (nondichotomic). The latter regime is achieved simply by scaling the gene number with Ω together with the rest of the protein and *mRNA* populations. We then compare the 2D probability distributions with the corresponding deterministic limit cycles (Eq. **5**). First, from the comparison with the deterministic limit cycles, we see that in both adiabatic and nonadiabatic cases, there are high-probability regions that lie outside the deterministic path, showing that the deterministic description does fully capture all the features of stochastic gene oscillators. Second, in single-gene cases, there are relatively large fluctuations along the radial directions, even in the near-deterministic limit Ω ∼ 10^{3}. Although the probability distribution is somewhat more localized around the deterministic cycle, in the adiabatic case, it still is not quite like the nondichotomic case, in which there is only a very narrow distribution around the deterministic cycle at high Ω > 10^{2} values. It is interesting to note that contrary to the nondichotomic case, both adiabatic and nonadiabatic oscillators never quite approach the deterministic limit, even as Ω → ∞. This indicates that gene oscillators do not have a precise deterministic analog.

By computing the autocorrelation function for various Ω values (Fig. 6), we see that the dephasing time reaches a plateau instead of diverging, as is the case with the well-investigated case of chemical oscillators with nondichotomic molecular noise. Thus, in the case of genetic oscillators, there is residual noise associated with the genetic degree of freedom, which is not completely eliminated by increasing protein and *mRNA* numbers. In the master equation language, this means that spectral gap is not a simple linear function of the system size Ω.

Besides making the oscillations more stochastic, the slowing of gene-state fluctuations also delays the negative feedback exerted by the *IκB*. The delay is reflected in the expansion of the limit cycle (Fig. S2 and Fig. 5) or lowering of the oscillation frequencies, as is seen from the systematic leftward shift of the power spectrum peak (29) (Fig. S4). Interestingly, the noise-induced expansion of the limit cycle is qualitatively captured by the phenomenological model (see the expression for ), which predicts higher oscillation amplitudes with more “noisy” oscillators. One expects the rate of dephasing to be a function of the oscillation cycling rate or the period. In other words, oscillators that individually complete their cycles faster given the same levels of noise also will dephase faster in the ensemble.

We can tune the cycling rate without significantly perturbing the noise levels by adjusting the irreversible rates that are part of the negative feedback loop of the network while keeping the binding/unbinding rates fixed. By doing this, one finds that both faster cycling rates and slower binding/unbinding rates do indeed lead to faster dephasing, and vice versa (Fig. 7). Lowering the gene binding/unbinding rates, we observed a steady increase in the noise level (Fig. S3), whereas changing translation or transcription rates had virtually no impact on the noise level. This confirms our reasoning that cycling rate modulates the dephasing time scale by propagating noise faster. Thus, the factors of phase noise and cycling rate are seen as the main contributors to dephasing, in which the former quantifies the rate of decorrelation but only the latter is needed to set the time scale of the oscillations.

It is expected that dephasing in vivo will be controlled by extrinsic cellular machinery via coupling of different oscillators. This not only should reduce the basal phase noise, but also may lead to a high degree of temporal organization of the underlying biochemical events. The role of coupling strength in the dephasing rates will be pursued in future study. Quantitative understanding of the principles of dephasing in gene oscillators should help in the interpretation and comparison of single-cell and population-based experiments and also may help in the design of synthetic gene oscillatory circuits.

## Conclusion

In the present work, two mathematical approaches are elaborated, each having its own unique advantages for understanding the dephasing phenomena in gene oscillators. Discrete-Markov models, although computationally demanding, provide a rigorous way to treat dephasing as a nonequilibrium relaxation from a synchronized state to the state of uniform phase distribution. Phenomenological models, although less rigorous, give valuable qualitative insights. However, care must be taken when using such models because the stochastic simulations on gene oscillators show that dichotomic noise from gene-state switching leads to significant deviations from the models that are rooted in the bifurcations presented in the deterministic equations of motion.

## Acknowledgments

We thank Dr. Jin Wang, Dr. Aleksandra Walczak, Dr. Masaki Sasai, and Dr. Bin Zhang for critical reading of the manuscript. We gratefully acknowledge financial support by the D.R. Bullard-Welch Chair at Rice University and PPG Grant P01 GM071862 from the National Institute of General Medical Sciences.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: pwolynes{at}rice.edu.

Author contributions: D.A.P. and P.G.W. designed research; D.A.P. performed research; D.A.P. analyzed data; and D.A.P. and P.G.W. wrote the paper.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1323433111/-/DCSupplemental.

Freely available online through the PNAS open access option.

## References

- ↵
- Winfree A

- ↵
- Schultz D,
- Wolynes PG,
- Ben Jacob E,
- Onuchic JN

- ↵
- ↵
- Goldbeter A

- ↵
- ↵
- ↵
- Hoffmann A,
- Levchenko A,
- Scott ML,
- Baltimore D

- ↵
- Ashall L,
- et al.

- ↵
- Nelson DE,
- et al.

- ↵
- Krishna S,
- Jensen MH,
- Sneppen K

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Wang J,
- Xu L,
- Wang E

- ↵
- ↵
- ↵
- ↵
- ↵
- Guckenheimer J,
- Holmes P

- ↵
- Stratonovich R,
- Silverman R

- ↵
- Baras F

- ↵
- ↵
- ↵
- ↵
- Bergqvist S,
- et al.

- ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Biological Sciences
- Systems Biology

- Physical Sciences
- Statistics