Chemical event chain model of coupled genetic oscillators

We introduce a stochastic model of coupled genetic oscillators in which chains of chemical events involved in gene regulation and expression are represented as sequences of Poisson processes. We characterize steady states by their frequency, their quality factor, and their synchrony by the oscillator cross correlation. The steady state is determined by coupling and exhibits stochastic transitions between different modes. The interplay of stochasticity and nonlinearity leads to isolated regions in parameter space in which the coupled system works best as a biological pacemaker. Key features of the stochastic oscillations can be captured by an effective model for phase oscillators that are coupled by signals with distributed delays.

Figure 1

Exemplary depiction of a simple genetic feedback system with a gene $x$ and a gene product $y$. (a) Simplified feedback scheme that can be mapped onto a dynamic system and (b) representation as a chemical chain involving intermediate chemical states $x_i$ and $y_i$ that can be mapped onto a master equation.

Figure 2

The zebrafish somitogenesis oscillator as an example for coupled genetic oscillations. Shown are two cells that act as autonomous oscillators and that are coupled through the Delta-Notch signaling pathway, an example for juxtacrine signaling [37, 49]. Coupling is bidirectional, that is, each cell acts as both a sender and a receiver. During embryonic development, a tissue comprising these cellular oscillators guides the segmentation of the elongating body axis [19].

Figure 3

(a) Schematics of the chemical event chain model of two coupled genetic oscillators as described by Eqs. (1) and (2). Boxes mark the initial and final products of a multistep process and broken lines the intermediate products. (b) Hill functions ${\psi }_{-}$ and ${\psi }_{+}$ as given by Eqs. (2) for different values of the exponent $h$.

Figure 4

Example time series of the cyclic final products ${\tilde{x}}_{1\tilde{n}}$ and ${\tilde{x}}_{2\tilde{n}}$ for different synchronization scenarios: (a) in-phase, (b) uncorrelated, and (c) antiphase. For all plots, parameters are given in Table 1 except for the transition rate $\lambda$, which is chosen such that the effective coupling delay $\tau$, given by Eq. (4), takes values $\tau =0.1T$ (a), $\tau =0.2T$ (b), and $\tau =0.5T$ (c), where $T=28$ is the period of the uncoupled oscillators.

Figure 5

Coupling delays determine the mode of synchrony and the collective frequency. [(a)–(c)] Density plots of the cross correlation $C$, Eq. (8). Blue colors indicate positive values of $C$ corresponding to in-phase correlations; red colors indicate negative values of $C$ corresponding to antiphase correlations. [(d)–(f)] Density plots of the collective frequency $\mathrm{\Omega }$. The parameters that are not varied are given in Table 1.

Figure 6

Coupling enhances precision. Density plots of the quality factor $Q$, Eq. (7), as a function of the scaled mean coupling delay $\tau /T$ and (a) the activation strength $\beta$, (b) the activation threshold $q$, and (c) the decay rate $\kappa$ of signaling molecules. The graded bar on top of each panel indicates the range of $Q/Q_0$ values where $Q_0=20.8$ is the quality of an uncoupled oscillator. The parameters that are not varied are given in Table 1.

Figure 7

Stochastic switching between antiphase and in-phase synchronized oscillations. (a) Time series of products ${\tilde{x}}_{1\tilde{n}}$ and ${\tilde{x}}_{2\tilde{n}}$ showing a switching event between antiphase and in-phase oscillations. Parameters are given in Table 1 except for $\beta =120$ and $\lambda =5/3$. (b) Windowed cross correlation $c$ of the time series in (a), see Eq. (9), showing the three distinct regimes. The width of the averaging window is $w=4T$. (c) Wavelet scalogram of one of the time series in (a). Bright regions indicate strong period components. The white dashed line serves as a guide for the eye.

Figure 8

Collective frequency $\mathrm{\Omega }$ of the phase oscillator system Eq. (11) (curves) and its asymptotic state for given constant initial phase differences $\mathrm{\Delta }\varphi$ (region plots) as a function of the scaled coupling delay $\tau /T$ (see Sec. 4b). Blue curves show the in-phase synchronized state, and red curves show the antiphase synchronized state. Solid curves show stable solutions, and dashed curves show unstable solutions. (a) System with a distribution of delay times for $n=10$ and different mean delays $\tau$ obtained by varying $\lambda$. The collective frequency is determined by Eq. (12). (b) Collective frequency for a discrete delay time $\tau$, determined by Eq. (14). In both plots, $\omega =0.224$ and $\varepsilon =\omega /4$.

Figure 9

Collective frequency as a function of the scaled coupling delay $\tau /T$ for different activation strengths $\beta$: [(a) and (e)] $\beta =30$, [(b) and (f)] $\beta =60$, [(c) and (g)] $\beta =90$, [(d) and (h)] $\beta =120$. [(a)–(d)] The density plots show logarithmic periodograms of oscillations in the chemical event chain model, where bright regions correspond to strong frequency components. [(e)–(h)] The curves (superimposed on the same density plots as in (a)–(d)] show the collective frequencies of the in-phase (black) and antiphase (white) solutions of the phase model, Eq. (12). Solid lines show stable solutions, and dashed lines show unstable solutions. The delay is given in multiples of the uncoupled period $T$, and the collective frequency $\mathrm{\Omega }$ is given in multiples of the uncoupled frequency $\omega =2\pi /T$. Solid lines indicate stable solutions, and dashed lines indicate unstable solutions. The parameters for the chemical event chain model are given in Table 1. The parameters for the phase model are adopted from the chemical event chain model with $\omega \approx 0.224$ estimated via Eq. (15) and $\varepsilon =r\beta$ with $r=1.33\times {10}^{-4}$.