Intrinsic noise, Delta-Notch signalling and delayed reactions promote sustained, coherent, synchronised oscillations in the presomitic mesoderm

TL;DR

This paper uses stochastic modeling to show how Delta-Notch signalling, intrinsic noise, and delays contribute to synchronized, sustained oscillations in the presomitic mesoderm, crucial for somite segmentation.

Contribution

It demonstrates that intrinsic noise and delays actively promote and enhance coherent oscillations, extending understanding beyond deterministic models.

Findings

Delta-Notch signalling synchronizes gene expression cycles.

Intrinsic noise amplifies and sustains oscillations.

Delays can cause anti-phase oscillations and stochastic wave propagation.

Abstract

Using a stochastic individual-based modelling approach, we examine the role that Delta-Notch signalling plays in the regulation of a robust and reliable somite segmentation clock. We find that not only can Delta-Notch signalling synchronise noisy cycles of gene expression in adjacent cells in the presomitic mesoderm (as is known), but it can also amplify and increase the coherence of these cycles. We examine some of the shortcomings of deterministic approaches to modelling these cycles and demonstrate how intrinsic noise can play an active role in promoting sustained oscillations, giving rise to noise-induced quasi-cycles. Finally, we explore how translational/transcriptional delays can result in the cycles in neighbouring cells oscillating in anti-phase and we study how this effect relates to the propagation of noise-induced stochastic waves.

Equations

$${\delta }_j^{\alpha }\left(t\right)=n_j^{\alpha }\left(t\right)-{\bar{n}}_j^{\alpha }\left(t\right).$$

$$P_j^{\alpha }\left(\omega \right)=⟨|{\hat{\delta }}_j^{\alpha }\left(\omega \right){|}^2⟩,$$

$$\tan \left({\varphi }_{j,j^{\prime }}^{\alpha ,{\alpha }^{\prime }}\left(\omega \right)\right)=\frac{\mathrm{Im}\left(⟨{\hat{\delta }}_j^{\alpha }\left(\omega \right){\hat{\delta }}_{j^{\prime }}^{{\alpha }^{\prime }\star }\left(\omega \right)⟩\right)}{\mathrm{Re}\left(⟨{\hat{\delta }}_j^{\alpha }\left(\omega \right){\hat{\delta }}_{j^{\prime }}^{{\alpha }^{\prime }\star }\left(\omega \right)⟩\right)}.$$

$$A_j^{\alpha }={\int }_0^{\infty }{P}_j^{\alpha }\left(\omega \right)d\omega .$$

$$C_j^{\alpha }=\frac{1}{A_j^{\alpha }}{\int }_{{\omega }_{\max }^{(j)}-\frac{\Delta \omega }{2}}^{{\omega }_{\max }^{(j)}+\frac{\Delta \omega }{2}}{P}_j^{\alpha }\left(\omega \right)d\omega .$$

$$P_k^{\alpha }\left(\omega \right)=⟨|{\hat{\tilde{\delta }}}_k^{\alpha }\left(\omega \right){|}^2⟩,$$

$$M_d$$

$$M_h$$

$$\varnothing$$

$$P_d$$

$$P_h$$

$$M_d$$

$$M_h$$

$$f_j^{\alpha }\left({\mathbf{n}}_t\right)$$

$$q_{j,t+\Delta }^{\alpha }-q_{j,t}^{\alpha }=\sum _{r,\tau }\frac{{\nu }_r^{\alpha }{k}_{j,r,\tau ,t}}{N}+\sum _{r,\tau }\frac{w_{\tau ,r}^{\alpha }{k}_{j,r,\tau ,t-\tau }}{N}.$$

$$P\left(\mathbf{n},t\right)=\frac{1}{N}\int \left(\prod _{j,\alpha }\prod _{t^{\prime }=0}^{t}d{q}_{j,t^{\prime }}^{\alpha }\right)\delta \left[\frac{1}{N}\mathbf{n}-\mathbf{q}\left(t\right)\right]\mathcal{P}\left(\{q_{j,t^{\prime }}^{\alpha }\}\right),$$

$$P\left(\mathbf{n},t\right)=\frac{1}{N}\sum _{\{k\}}\int \left(\prod _{j,\alpha }\prod _{t^{\prime }=0}^{t}d{q}_{j,t^{\prime }}^{\alpha }\right)\delta \left[\frac{1}{N}\mathbf{n}-\mathbf{q}\left(t\right)\right]\prod _{j,\alpha ,\tau ,t^{\prime }}\delta \left(q_{j,t^{\prime }+\Delta }^{\alpha }-q_{j,t^{\prime }}^{\alpha }-\sum _{r,\tau }\frac{{\nu }_r^{\alpha }{k}_{j,r,\tau ,t^{\prime }}}{N}-\sum _{r,\tau }\frac{w_{\tau ,r}^{\alpha }{k}_{j,r,\tau ,t^{\prime }-\tau }}{N}\right)P\left(\{k\}\right)$$

$$=\frac{1}{N}\sum _{\{k\}}\int \left(\prod _{j,\alpha }\prod _{t^{\prime }=0}^t\frac{d{q}_{j,t^{\prime }}^{\alpha }d{p}_{j,t^{\prime }}^{\alpha }}{2\pi }\right)\delta \left[\frac{1}{N}\mathbf{n}-\mathbf{q}\left(t\right)\right]e^{i{\sum }_{j,\alpha ,t^{\prime }}{p}_{j,t^{\prime }}^{\alpha }\left(q_{j,t^{\prime }+\Delta }^{\alpha }-q_{j,t^{\prime }}^{\alpha }-{\sum }_{r,\tau }\frac{{\nu }_r^{\alpha }{k}_{j,r,\tau ,t^{\prime }}}{N}-{\sum }_{r,\tau }\frac{w_{\tau ,r}^{\alpha }{k}_{j,r,\tau ,t^{\prime }-\tau }}{N}\right)}P\left(\{k\}\right),$$

$$P\left(\{k\}\right)=\prod _{j,r,\tau ,t}P\left(k_{j,r,\tau ,t}\right).$$

$$P\left(k_{j,r,\tau ,t}\right)=\frac{{\left({\Delta }^2{W}_{j,r,\tau ,t}\right)}^{k_{j,r,\tau ,t}}}{k_{j,r,\tau ,t}!}{e}^{-{\Delta }^2{W}_{j,r,\tau ,t}}.$$

$$\sum _{k_{j,r,\tau ,t}}\frac{{\left({\Delta }^2{W}_{j,r,\tau ,t}\right)}^{k_{j,r,\tau ,t}}}{k_{j,r,\tau ,t}!}{e}^{-{\Delta }^2{W}_{j,r,\tau ,t}}{e}^{-i{\sum }_{\alpha }{p}_{j,t}^{\alpha }\frac{{\nu }_r^{\alpha }{k}_{j,r,\tau ,t}}{N}-i{\sum }_{\alpha }{p}_{j,t+\tau }^{\alpha }\frac{w_{\tau ,r}^{\alpha }{k}_{j,r,\tau ,t}}{N}}$$

$$=\exp \left\{-{\Delta }^2{W}_{j,r,\tau ,t}\left[1-e^{-i{\sum }_{\alpha }\frac{{\nu }_r^{\alpha }{p}_{j,t}^{\alpha }+w_{\tau ,r}^{\alpha }{p}_{j,t+\tau }^{\alpha }}{N}}\right]\right\}.$$

$$P\left(\mathbf{n},t\right)=\frac{1}{N}\int \left(\prod _{j,\alpha }\prod _{t^{\prime }=0}^t\frac{d{q}_{j,t^{\prime }}^{\alpha }d{p}_{j,t^{\prime }}^{\alpha }}{2\pi }\right)\delta \left[\frac{1}{N}\mathbf{n}-\mathbf{q}\left(t\right)\right]e^{i{\sum }_{j,\alpha ,t^{\prime }}{p}_{j,t^{\prime }}^{\alpha }\left(q_{j,t^{\prime }+\Delta }^{\alpha }-q_{j,t^{\prime }}^{\alpha }\right)}$$

$$\times \exp \left\{-\sum _{j,r,\tau ,t^{\prime }}{\Delta }^2{W}_{j,r,\tau ,t^{\prime }}\left[1-e^{-i{\sum }_{\alpha }\frac{{\nu }_r^{\alpha }{p}_{j,t^{\prime }}^{\alpha }+w_{\tau ,r}^{\alpha }{p}_{j,t^{\prime }+\tau }^{\alpha }}{N}}\right]\right\}.$$

$$P\left(\mathbf{n},t\right)=\frac{1}{N}\int \left(\prod _{j,\alpha }\prod _{t^{\prime }=0}^t\frac{d{q}_{j,t^{\prime }}^{\alpha }d{p}_{j,t^{\prime }}^{\alpha }}{2\pi }\right)\delta \left[\frac{1}{N}\mathbf{n}-\mathbf{q}\left(t\right)\right]e^{i{\sum }_{j,\alpha ,t^{\prime }}{p}_{j,t^{\prime }}^{\alpha }\left(q_{j,t^{\prime }+\Delta }^{\alpha }-q_{j,t^{\prime }}^{\alpha }\right)}$$

$$\displaystyle\times\exp\Bigg{\{}-\sum_{j,r,\tau,t^{\prime}}\Delta^{2}W_{j,r,\tau,t^{\prime}}\Bigg{[}\frac{i}{N}\sum_{\alpha}\left(\nu_{r}^{\alpha}p^{\alpha}_{j,t^{\prime}}+w_{\tau,r}^{\alpha}p^{\alpha}_{j,t^{\prime}+\tau}\right)$$

$$\displaystyle+\frac{1}{2N^{2}}\sum_{\alpha,\alpha^{\prime},j^{\prime},t^{\prime},\tau^{\prime}}\delta_{\tau,\tau^{\prime}}\delta_{j,j^{\prime}}\delta_{t^{\prime},t^{\prime\prime}}\left(\nu_{r}^{\alpha}p^{\alpha}_{j,t^{\prime}}+w_{\tau,r}^{\alpha}p^{\alpha}_{j,t^{\prime}+\tau}\right)\left(\nu_{r}^{\alpha^{\prime}}p^{\alpha^{\prime}}_{j^{\prime},t^{\prime\prime}}+w_{\tau,r}^{\alpha^{\prime}}p^{\alpha^{\prime}}_{j^{\prime},t^{\prime\prime}+\tau^{\prime}}\right)\Bigg{]}\Bigg{\}},$$

$$\frac{d{n}_{j,t}^{\alpha }}{dt}$$

$$⟨{\xi }_{j,t}^{\alpha }{\xi }_{j^{\prime },t^{\prime }}^{{\alpha }^{\prime }}⟩$$

$$P\left(\mathbf{n},t\right)=\frac{1}{N}\int \left(\prod _{j,\alpha }\prod _{t^{\prime }=0}^t\frac{d{q}_{j,t^{\prime }}^{\alpha }d{p}_{j,t^{\prime }}^{\alpha }}{2\pi }\right)\delta \left[\frac{1}{N}\mathbf{n}-\mathbf{q}\left(t\right)\right]e^{i{\sum }_{j,\alpha ,t^{\prime }}{p}_{j,t^{\prime }}^{\alpha }\left(q_{j,t^{\prime }+\Delta }^{\alpha }-q_{j,t^{\prime }}^{\alpha }\right)}$$

$$\displaystyle\times\exp\Bigg{[}-\frac{i\Delta}{N}\sum_{\alpha,j,t^{\prime}}p^{\alpha}_{j,t^{\prime}}F^{\alpha}_{j,t^{\prime}}-\frac{\Delta^{2}}{2N^{2}}\sum_{\alpha,\alpha^{\prime},j,j^{\prime},t^{\prime},t^{\prime\prime}}p^{\alpha}_{j,t^{\prime}}p^{\alpha^{\prime}}_{j^{\prime},t^{\prime\prime}}\Sigma^{\alpha\alpha^{\prime}}_{j,j^{\prime},t^{\prime},t^{\prime\prime}}\Bigg{]}.$$

$$\frac{d{n}_{j,t}^{\alpha }}{dt}$$

$$⟨{\xi }_{j,t}^{\alpha }{\xi }_{j^{\prime },t^{\prime }}^{{\alpha }^{\prime }}⟩$$

$$+$$

$$\frac{d{n}_{j,t}^{\alpha }}{dt}$$

$$+\sum _r{\int }_0^{\infty }{W}_{j,r,t-\tau }{K}_r\left(\tau \right)w_{r,\tau }^{\alpha }d\tau +{\xi }_{j,t}^{\alpha }.$$

$$\frac{d{\delta }_{j,t}^{\alpha }}{dt}$$

Full text

I Introduction

In developing vertebrates and cephalochordates, as the embryo forms and extends pairs of blocks of mesodermal progenitor cells assemble, bilaterally flanking the notochord Pourquié (2003). These blocks, termed somites, eventually go on to form vertebrae and ribs after further cellular differentiation. The somites are constructed pair-by-pair, anterior to posterior, in a rythmic and sequential manner as the tailbud extends away from the rostral end of the embryo. They are formed from cells originating in the presomitic mesoderm (PSM). Such cells are produced continually by the tailbud as the abdomen elongates Oates et al. (2012); Saga and Takeda (2001).

The process of somite segmentation has been of interest to experimentalists and theorists working in the field of developmental biology for some decades; it provides a fascinating case study where one can directly examine the link between microscopic gene regulatory systems operating in individual cells and macroscopic developmental processes. The prevailing theoretical framework for understanding the process was put forward by Cooke and Zeeman in 1976 Cooke and Zeeman (1976) and is termed the ‘clock-wavefront’ model. This model proposes that the cells in the presomitic mesoderm each possess an internal cyclic ‘clock’ which is synchronised between the cells. Additionally, a wavefront propagates through the PSM as the embryo grows. As the wavefront encounters cells, it interacts with them differently depending on the current state of each internal cellular clock. This interaction causes the cells to change their adhesive and migratory properties. The temporal periodicity of the cell cycles is thus converted into the spatial periodicity of the somites.

Considerable experimental and theoretical effort has been expended in order to identify the genetic oscillators that constitute the putative somite segmentation ‘clock’ and a good amount of progress has been made. In certain model organisms, such as the mouse and the zebrafish, so-called ‘knockdown’/‘knockout’ experiments have identified genes which when mutated give rise to defects in the formation of somites and, consequently, the vertebrae Bessho et al. (2001); Henry et al. (2002); Oates and Ho (2002); Sieger et al. (2006); Conlon et al. (1995); de Angelis et al. (1997); Kusumi et al. (1998). Gradients of FGF (fibroblast growth factor) or Wnt protein, which are produced in the tailbud, are thought to constitute the moving wavefront; transient loss or increase in these substances can alter the local somite length Aulehla et al. (2003); Dubrulle et al. (2001). Genes such as hes in the mouse Jouve et al. (2000) and her in the zebrafish Holley et al. (2000) are thought to be the primary cyclic genes which act as clocks. These genes are both targets of the Notch signalling pathway. It has also been shown in experiments that Delta-Notch signalling is a vital component in synchronising oscillations Jouve et al. (2000); Jiang et al. (2000); Özbudak and Lewis (2008); Delaune et al. (2012).

In order for the oscillations in the expression of the hes/her genes to constitute a viable segmentation clock for the clock-wavefront model, the oscillations must satisfy several criteria: (1) The oscillations in gene expression must have the same frequency in adjacent cells. (2) The oscillations in adjacent cells must be in phase. (3) The cellular oscillations must be coherent (there must be a clear dominant frequency). (4) The oscillations must have sizeable enough an amplitude so as not be indistinguishable from background ‘noise’. Mathematical models of the gene regulatory system have shown that Delta-Notch signalling can indeed synchronise (align the frequencies) of oscillations in neighbouring cells with intrinsically differing cellular clocks Lewis (2003). That is, it has been shown that Delta-Notch signalling is responsible for satisfying condition (1), but relatively little discussion has been dedicated to the latter 3 conditions (phase, coherence and amplitude).

So that one might analyse the degree to which the criteria above are satisfied, one must take into account stochastic (random) effects in the system, especially with regards to point (3). The gene regulatory systems in question are inherently noisy in nature Özbudak et al. (2002); Austin et al. (2006); Shimojo et al. (2016); Phillips et al. (2016, 2017); Hirata et al. (2002); Shimojo et al. (2008). This is due in part to the stochastic nature of the production/decay events of individual proteins and/or mRNA molecules and the fact that there are finite numbers of these molecules in any one cell. Noise of this origin is termed intrinsic in the literature Elowitz et al. (2002); Swain et al. (2002). On the other hand, gene expression is also influenced by the concentrations, locations and states of molecules such as regulatory proteins or polymerases which can affect the global activity in a single cell but can vary between cells. Noise arising from fluctuations in the properties of such molecules is referred to as extrinsic Elowitz et al. (2002); Swain et al. (2002).

The role of noise has largely been disregarded in previous theoretical work on the somite segmentation clock Jensen et al. (2010); Goldbeter and Pourquié (2008); Momiji and Monk (2009), or has often been treated only as an external influence rather than as an aspect intrinsic to the translation and transcription processes Cinquin (2003); Morelli et al. (2009); Horikawa et al. (2006) (an exception can be found in Ay et al. (2014), where simulations involving intrinsic noise were performed). For example, some works have considered the binding and unbinding of repressor protein to the DNA binding site as a stochastic process Lewis (2003); Baker et al. (2008) so that the rates of transcription are themselves stochastic variables. This source of noise is taken into account in the context of deterministic evolution equations - the intrinsic stochasticity of the transcription and translation events is not accounted for. In this paper however, we study the effect of intrinsic stochasticity by treating the production and the decay of individual molecules as random processes, following Özbudak et al. (2002); Galla (2009); Barrio et al. (2006); Schlicht and Winkler (2008); Turner et al. (2004); Thattai and Van Oudenaarden (2001). These events may be subject to delays arising from the finite time taken for the translation/transcription processes.

Using an individual-based mathematical model capturing the intrinsic noise in the system, we demonstrate that, perhaps counter to intuition, intrinsic stochasticity can be a proactive force in promoting cellular oscillations. We study how these noisy oscillations in neighbouring cells are affected by different levels of the strength of Delta-Notch signalling. We are able to show that, under certain conditions, Delta-Notch signalling not only acts to align the frequencies of oscillations in neighbouring cells; it can also reduce the phase lag, reduce the range of dominant oscillatory frequencies (i.e., it can make the oscillations more coherent) and it can increase the amplitude of oscillations (also noted in Ay et al. (2014); Özbudak and Lewis (2008)). The combination of these effects indicates that Delta-Notch signalling can contribute to satisfying points (1)-(4) above. We also discuss circumstances under which pairs of cells may oscillate out-of-phase, despite Delta-Notch coupling. We explore how this is related to waves and/or oscillating chequerboard patterns of gene expression in extended chains of cells.

II Methods

II.1 Model definition

Oscillations in the expression of genes (or ‘pulsing dynamics’) is a well-documented phenomenon responsible for many cell functions and broader biological processes Levine et al. (2013); Goldbeter (1997). Such oscillations in the expression of certain genes are thought to constitute the biological ‘clock’ in the clock-wavefront model Cooke and Zeeman (1976) of somite segmentation Saga and Takeda (2001); Oates et al. (2012). The genes in question are known to be affected by Notch signalling. For the purposes of our theoretical treatment it is not necessary to consider the full complexity of the Notch signalling pathway Bray (2006); Andersson et al. (2011) or even the full network of interacting genes involved with the somite segmentation process Pourquié (2011); one can use a simplified model to highlight the salient features and analyse their causes.

So-called ‘knockdown’/‘knockout’ experiments Saga and Takeda (2001) suggest that the most relevant genes for the regulation of the ‘clock’ are delta (or its homologues) and hes in mice, and her in zebrafish Oates et al. (2012). It has been shown previously Momiji and Monk (2009) that a two-gene model involving only hes/her and delta is sufficient for the emergence of cycles. We therefore also adopt a reduced two-gene model, as this will be sufficient to highlight the effects with which we are concerned. The reduced system is depicted schematically in Fig. 1 and is discussed in more detail in the Supplement (Section S1). For now, we consider a system of two coupled cells as a simple example. We generalise the approach to systems of greater numbers of cells in Section III.4.

Using $\mathbf{n}\left(t\right)$ to denote the set of all protein and mRNA numbers in all cells at time $t$, the reduced gene regulatory system that we consider can be summarised as follows: hes/her mRNA is transcribed at a rate $f^{h}\left[\mathbf{n}\left(t-\tau^{(m)h}\right)\right]$. The rate $f^{h}\left[\mathbf{n}\left(t-\tau^{(m)h}\right)\right]$ takes the form of a sum of Hill functions, which reflects the fact that hes/her mRNA production is inhibited by local Hes/Her protein and activated by Delta protein in adjacent cells. The precise form of this non-linear function is given in the Supplement (Section S1). Every hes/her mRNA molecule is translated into protein at a constant rate $a_{h}$. Because the transcription and translation processes take an amount of time $\tau^{(m)h}$ and $\tau^{(p)h}$ to complete respectively, the present rate of production of mRNA/protein is dependent upon protein/mRNA concentrations in the past (respectively). Both hes/her mRNA and protein molecules degrade (become inert) at constant per capita rates $c_{h}$ and $b_{h}$ respectively. In a similar way, delta mRNA is transcribed at a rate $f^{d}\left[\mathbf{n}\left(t-\tau^{(m)d}\right)\right]$ and decays at a constant per capita rate $c_{d}$. Delta protein molecules are produced and decay at the per capita rates $a_{d}$ and $b_{d}$ respectively. Production of delta mRNA and protein are delayed by the times $\tau^{(m)d}$ and $\tau^{(p)d}$ to respectively. Finally, production of delta mRNA is inhibited by local Hes/Her protein concentration.

There are three vital aspects to the processes in this setup: (1) The model is individual-based – it does not treat protein/mRNA concentrations as continuous quantities (an approximation only valid when population numbers are large). The production and degradation of proteins and mRNA are inherently stochastic (random) processes due to the finite numbers of proteins and mRNA Phillips et al. (2016, 2017) in each cell; this gives rise to noisy dynamics Turner et al. (2004). (2) There is a time-delay between the activation of the production of one unit of mRNA/protein and the completion of the production process. As a result, the rates of production of mRNA/protein at a given time are dependent on the state of the system in the past. In the language of stochastic processes, the dynamics are non-Markovian (they have memory) Gardiner (2009). It has been established that time-delays such as these can encourage the emergence of temporal oscillations Momiji and Monk (2009); Lewis (2003); Smith (2011). (3) Due to Delta-Notch signalling, the rate of production of hes/her mRNA in one cell is dependent on the concentration of Delta protein in the neighbouring cells. In this sense, there is a non-locality to the reaction rates.

The combination of these three aspects of the dynamics leads to a unique challenge with respect to theoretical modelling. However, we demonstrate in the Supplement that one can approximate the full individual-based dynamics of the system with a set of stochastic differential equations (SDEs); these are given in the Supplement (Section S2 C). They take a similar form to the deterministic (noiseless) equations given in Lewis (2003) but include additional Gaussian noise terms which take into account the intrinsic stochasticity of the system. We emphasize that the properties of this noise are calculated so as to agree with individual-based simulations of the system; the noise is not added in an ad hoc fashion. The tools we use to quantify the phenomena induced in the system by noise are discussed in the next section.

II.2 Analysis of stochastic behaviour

In this work, we will be concerned primarily with the theoretical analysis of noise-induced cycles of gene expression. These are oscillations which occur in the full stochastic individual-based model but which are missing in the noiseless deterministic system.

The power spectrum of fluctuations about the deterministic trajectory will be the main quantitative tool that we use to analyse the noise-induced phenomena in the gene regulatory model described in Section II.1 (and elaborated upon in the Supplement Section S1). We denote the number of particles of type $\alpha$ in cell $j$ at time $t$ by $n^{\alpha}_{j}(t)$. The type of particle indicated by the index $\alpha$ may be mRNA molecules or proteins. The dynamics of the quantities $n^{\alpha}_{j}(t)$ are approximated by the system of stochastic differential equations (SDEs) in Eq. (S25) (in the Supplement). Further, we write $\bar{n}^{\alpha}_{j}\left(t\right)$ for the numbers of particles predicted by the corresponding deterministic model [Eq. (S25), with the noise terms $\xi^{\alpha}_{j}\left(t\right)$ set to zero]. Thus, we define the fluctuations about the deterministic trajectory as

[TABLE]

The power spectrum of fluctuations is then defined via the temporal Fourier transform as

[TABLE]

where the Fourier transform is given by $\hat{g}\left(\omega\right)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}e^{i\omega t}g\left(t\right)dt$; the angular brackets denote the ensemble average over the set of all possible stochastic time courses of the system. Roughly speaking, the power spectrum of fluctuations decomposes a time series into its composite frequencies and quantifies the statistical contribution of a particular frequency to the series. A large, narrow, unique peak in the power spectrum indicates that the frequency at which the peak is located is the dominant frequency of the time series; a peak of this type centred on a non-zero frequency therefore characterises periodicity. If the deterministic trajectory $\bar{n}^{\alpha}_{j}\left(t\right)$ is non-oscillatory, then such a peak in the spectrum of the stochastic model indicates noise-induced oscillations. We note that in our analysis the power spectrum is always evaluated in the steady state, i.e. when all transient effects have decayed sufficiently so as to be negligible.

Using a mathematical approach (see Supplement Sections S1, S2 and S3), we are able to predict the power spectrum of the fluctuations when the deterministic trajectory has reached a fixed point (i.e., when $\bar{n}^{\alpha}_{j}\left(t\right)=\bar{n}^{\alpha}_{j}$ is constant at long times $t$). This theoretical approach gives us a way to identify, without performing time-consuming simulations, what the dominant frequency of noise-induced oscillations is and to what extent the other frequencies contribute. This analysis is only valid for noise-induced cycles, i.e. when the oscillations of the deterministic equations are transient. For the parameter regimes where there are persistent cycles in the deterministic equations, we perform a linear stability analysis in order to obtain quantities such as periods of oscillation and phase lags (see Supplement Section S3).

It is also possible to quantify the phase lag between two sets of noise-induced cycles in coupled cells using our theoretical approach. Following Challenger and McKane (2013); Rozhnova et al. (2012), we define the phase lag $\phi^{\alpha,\alpha^{\prime}}_{j,j^{\prime}}\left(\omega\right)$ associated with a particular frequency $\omega$ between species $\alpha$ in cell $j$ and species $\alpha^{\prime}$ in cell $j^{\prime}$ as

[TABLE]

Phases differing by integer multiples of $2\pi$ are degenerate therefore, in this paper, we define the phase lag to be in the range $[-\pi,\pi)$. Notably, the phase lag as defined in Eq. (3) is dependent on $\omega$. As mentioned above, a time series can be thought of as being comprised of a sum of cycles with different frequencies $\omega$. The quantity $\phi^{\alpha,\alpha^{\prime}}_{j,j^{\prime}}\left(\omega\right)$ is the phase lag between the constituent cycles of frequency $\omega$ in cells $j$ and $j^{\prime}$. In our analysis, we may refer to the phase lag of a cell $j$ (with respect to another cell), which we define as the phase lag at the peak frequency of the power spectrum of the cell in question, $\omega^{(j)}_{\textrm{max}}$.

Furthermore, following Alonso et al. (2006), we define the total amplification of fluctuations for particles of type $\alpha$ in cell $j$

[TABLE]

This quantity is proportional to the time-averaged squared displacement of the dynamics from the fixed point, i.e., to the variance of the stochastic time series.

Finally, again following Alonso et al. (2006), we also define the coherence as the proportion of the power spectrum within a fixed range $\Delta\omega$ of the peak

[TABLE]

The coherence $C^{\alpha}_{j}$ quantifies how sharply peaked the power spectrum is - i.e., how narrow the band of dominant frequencies is. It has a maximum value of $1$ and a minimum value of [math]. The choice of $\Delta\omega$ is largely immaterial provided $\Delta\omega$ is small compared to the peak frequency.

III Results

III.1 Individual-based models capture noise-driven effects which are missed by deterministic models

In the systems we are considering, individual cells contain of the order of 10-100 mRNA molecules and around 1000 proteins of any one type Phillips et al. (2016, 2017) (see also Lewis (2003); Keskin et al. (2018)). As such, the dynamics are inherently noisy. This type of noisy dynamics has been observed in experiments monitoring gene expression Austin et al. (2006); Shimojo et al. (2016); Phillips et al. (2016, 2017); Hirata et al. (2002); Shimojo et al. (2008). The expression of these genes cannot be fully described by the regular, smooth oscillation obtained from integrating deterministic sets of ODEs (as can be seen from Figs. 2 and 3). Instead, a stochastic individual-based model is better suited to qualitatively reproduce the results of experiment.

In previous theoretical studies of the somite segmentation clock, noise has mostly been treated as external to the dynamics Lewis (2003); Morelli et al. (2009); Horikawa et al. (2006) or has not been considered at all Momiji and Monk (2009); Kageyama et al. (2012). An exception to this is Ay et al. (2014), in which individual-based simulations were carried out (the consequences of the inclusion of intrinsic noise that we discuss here were not the focus of Ay et al. (2014) however). The noise that we use in the present work is rigorously derived as an intrinsic quality of the stochastic, individual-based dynamics themselves. As such, the results of our analysis agree with fully individual-based simulations of the system (as is demonstrated in Fig. 4).

It has previously been observed Lewis (2003) that noisy external driving can give rise to sustained oscillations in the somite segmentation clock. We show that similar noise-induced oscillations can also be produced by the stochastic nature of the transcription and translation processes in the segmentation clock themselves. That is, the inclusion of intrinsic noise in the theoretical modelling gives rise to sustained noise-induced oscillations which a purely deterministic model, or a more ad hoc approach to noise-inclusion, might miss. As is shown in Figs. 2 and 3, for sets of parameters which are biologically reasonable (see Lewis (2003)), one may observe the noiseless model tend towards a stationary fixed–point, only exhibiting transient oscillations which eventually decay. In the corresponding individual-based model however, the noise repeatedly ‘kicks’ the system away from the fixed point. As a consequence, the oscillations which were transient in the deterministic model are sustained by the noise. One thus observes persistent noisy oscillations for the same parameter set.

The oscillations in gene expression observed experimentally may very well be noise-induced cycles of this type. The emergence of such ‘quasi-cycles’ is a well-documented phenomenon which has been previously studied in the context of gene-regulatory models Galla (2009); Barrio et al. (2006) as well as in ecological systems Lugo and McKane (2008); McKane and Newman (2005) and in epidemics Alonso et al. (2006); Rozhnova et al. (2012); Black and McKane (2010). That these are indeed cycles with a periodic nature and not just random white noise is demonstrated by the power spectra of fluctuations (see Fig. 4)– this matter is discussed further in Sections II.2 and III.2.1.

Our theoretical approach to analysing noise-induced cycles (which is similar to that found in Brett and Galla (2013, 2014); Baron and Galla (2019) and detailed in the Supplement) allows us to study the amplification, synchronisation and coherence of these cycles, as discussed in the following sections.

III.2 Delta-Notch signalling mitigates inhomogeneity and promotes a robust and reliable segmentation clock in noisy oscillators

Having introduced the concept of noise-induced cycles and the mathematical tools that we will use to analyse them, we now turn our attention to the effect that increasing the Delta-Notch signalling strength has on these noisy oscillations.

It has been observed experimentally Jouve et al. (2000) that mutations in the delta gene give rise to defects in the formation of somites. This has been attributed to a decreased coupling between the cells arising from the mutation which, due to slight inhomogeneities between cells and the stochastic nature of the cellular cycles, leads to the genetic oscillations in neighbouring cells becoming asynchronous Jiang et al. (2000); Delaune et al. (2012). Furthermore, it has been shown that encumbered Delta-Notch signalling (i.e. reduced signalling strength) can give rise to greater disparities between the oscillations in cells which would be synchronised if signalling were not impaired Keskin et al. (2018). In this section, we reproduce and study this effect with our model and theoretical approach (presented in the Supplemental material in Section S2), thus verifying the necessity of Delta-Notch signalling for the somite segmentation clock. A justification for our mathematical definition of ‘coupling strength’ is given in Supplement Section S1.

Consider a two-cell system in which each cell has slightly different internal parameters (e.g. transcriptional delay time) such that the typical cycle time varies between the cells when they are uncoupled. We evaluate, using our theoretical approach, the response of the peak frequencies, inter-cell phase lag, amplification and coherence of genetic oscillations in this inhomogeneous two-cell system for various degrees of Delta-Notch coupling strength and thus show that the quality of the oscillations (and therefore the segmentation clock) can improve when Delta-Notch signalling is enhanced.

We use two different theoretical approaches in the following sections, each of which is valid for different, but complementary, parameter sets. (1) In the regime where the deterministic (noiseless) equations approach a fixed point, we evaluate the power spectrum of fluctuations of the emergent noise-induced cycles using the so-called linear-noise approximation (discussed in more detail Section S2 of the Supplemental Material). This analysis relies upon the deterministic dynamics tending towards a fixed point, and approximates the stochastic equations as linear in the vicinity of this fixed point. The accuracy of the approximation is tested against individual-based simulations of the full non-linear stochastic model in Fig. 4. (2) When the fixed point of the deterministic system becomes unstable, we use a deterministic linear stability analysis (LSA) to find the dominant oscillatory frequency of the cycles and the inter-cell phase lag. The LSA provides no way of finding the amplitude or the coherence of the cycles however – it is only possible for us to evaluate these when method (1) is valid.

It has been suggested that a sufficiently strong non-linearity (the so-called ‘cooperativity’) is required for regular sustained cycles Barrio et al. (2006). Noting that our linear analysis agrees with simulation results (see Fig. 4), we stress that the precise form of the non-linearity is not directly important for the emergence of noise-induced cycles. Their properties are well captured by the linearised dynamics. Non-linearities will however affect the location and nature of deterministic fixed points, and the coefficients in the linearised equations near these fixed points.

III.2.1 Delta-notch signalling synchronises noisy genetic oscillators and reduces their phase difference

Firstly, we find that Delta-Notch signalling can have the effect of synchronising the dominant oscillatory frequencies of two cells with differing internal parameters. This is demonstrated in Fig. 4, which depicts the power spectra of the stochastic fluctuations in two such cells for various degrees of coupling strength. One observes that as the coupling strength is increased, the peaks for either cell, which are separated when there is no coupling, are both drawn towards a common frequency, indicating synchronisation. The degree of synchronisation varies smoothly with the variation of the coupling strength, as shown in Fig. 5(a); in this figure, the dependence of the dominant frequencies on coupling strength is shown in more detail. The common frequency which is converged upon at large strengths of the Delta-Notch coupling agrees with that predicted by linear stability analysis (LSA) (see Supplement Section S3 and Section III.2.2 for further details).

Secondly, we find that the peaks of the power spectra in either cell converge to a frequency which reduces the phase lag between the oscillations in the two cells (see Fig. 4). So, not only can Delta-Notch signalling act to align the frequency of oscillations in neighbouring cells, it can also encourage the oscillators in either cell to be more aligned in phase. The smooth decrease of the phase lag with increasing coupling strength is shown in Fig. 5(b). In a similar way to the peak frequency, the phase lag between cells agrees with that predicted by LSA when the coupling is large.

Both of these factors, a shared oscillatory frequency and a minimal phase lag, are important for the proper functioning of a cellular clock. The changes in the peak frequencies and the phase lag that result from an increase in coupling strength correspond to quite a noticeable change in the quality of the oscillations themselves. Figs. 2 and 3 are evaluated for the same sets of parameters as Figs. 4(a) and 4(c) respectively. One parameter set is without coupling between the cells (Fig. 2) and one is with cell-to-cell coupling (Fig. 3). In the former case, the oscillations in either cell are somewhat aperiodic and there is no noticeable synchronisation between the cells. However, in Fig. 3, the highs and lows of the cycles in either cell are more inclined to align – this is associated with reduced phase lag and synchronised peak frequencies.

III.2.2 Delta-Notch coupling increases the amplitude and coherence of noisy oscillations

We asserted previously that an important characteristic of an effective cellular clock is a well-defined time-period – if many frequencies contribute significantly to the oscillations, then it is more difficult to identify an overall phase for the clock. We also asserted the necessity of the cycles to have a significant amplitude. Both of these factors contribute to the clarity of the ‘signal’ of the oscillations that constitute the cellular clock. We note that amplification and coherence are properties of the cycles in individual cells whereas synchronisation and phase lag are comparative measures of the oscillations in different cells. Despite this, amplification and coherence are indeed affected by Delta-Notch coupling too.

We find that as the Delta-Notch coupling strength is increased, the amplitude of the oscillations in both cells first decreases slightly then increases, as shown in Fig. 5(c). A similar result was also found in previous experimental and theoretical works Ay et al. (2014); Özbudak and Lewis (2008). However, it is important to note the following caveat on the results presented in Fig. 5(c):

The blue and red lines shown in Fig. 5 were produced using the theoretical approach based on the linear-noise approximation (see Sec. S2 of the Supplemental material). As such, we observe that the calculated amplitude of the oscillations diverges when the coupling strength is increased sufficiently. This singularity is a consequence of our theoretical approximation; it corresponds to the onset of an instability for the fixed point of the deterministic (noiseless) equations and the emergence of a limit cycle Boland et al. (2009, 2008). The nonlinearities in the model then become more relevant and curtail the amplitude of the oscillations; this is not captured by the linear theoretical approach. The point at which this instability occurs is predicted by deterministic LSA (see Supplement Section S3) and is indicated by dashed vertical lines in Fig. 5. It is at this point that the stochastic theory, which is valid only when the deterministic system approaches a fixed point, becomes inaccurate. Instead, the LSA becomes the more accurate tool for identifying the peak frequency of oscillation and the inter-cell phase lag. The two analytical methods complement each other in this sense – we can use both approaches to continuously analyse the dominant frequency and phase-lag over the onset of the deterministic instability. Unfortunately, the LSA provides no means of calculating the amplification or the coherence of the cycles – we are only able to accurately predict these quantities when the deterministic dynamics approaches a fixed point (i.e. before the onset of the deterministic instability).

With this caveat in mind, one can nevertheless see from Fig. 4 that the linear stochastic theory still agrees with simulation results close to this transition. That is, it remains accurate over a wide enough range of coupling strengths to faithfully capture an increase in the amplification with coupling strength, which occurs before the onset of the instability.

We find also that as the coupling strength is increased, the power spectrum of fluctuations becomes sharply peaked (as can be seen in Fig. 4(c)) at a characteristic frequency, i.e., the cycles become more coherent (see Fig. 5(d)). The location of this peak in the power spectrum corresponds to the frequency to which the cycles in the two cells converge (as discussed in the previous section).

We conceptualise this increase in amplification and coherence as a consequence of a kind of ‘resonant amplification’. Because of the communication between the cells, one can think of the state of one cell as influencing or ‘forcing’ the oscillations in the neighbouring cells. As the inter-cell coupling strength is increased, the frequencies are aligned and the phase lag between them is reduced, the cycles begin to constructively interfere at a characteristic frequency.

Interestingly, as the coupling strength is increased from zero, the amplitude of the oscillations in either cell initially decreases (Fig. 5(c)), as does the peakedness of the power spectrum (Fig. 5(d)). We attribute this to the fact that, for low coupling, the oscillations in either cell are not adequately synchronised for their interference to be constructive. So as the coupling strength is increased initially, the ‘interference’ between the two cells has a destructive effect. It is only when the phase lag is reduced and the frequencies are aligned sufficiently (as a result of a further increase of the coupling) that the collective amplitude of oscillations increases.

The effects of increased amplification and coherence are evident in the time-series of the noisy cellular oscillations in the two-cell system shown in Figs. 2 and 3 (which are evaluated with and without inter-cellular coupling respectively). In Fig. 2, there is no easily discernible periodic nature to the cycles in either cell. This is in contrast to the cycles in Fig. 3 where the highs and lows have more consistent temporal separations. This greater clarity in the oscillatory frequency is associated with the the increase in the sharpness of the peaks of the power spectra between Figs. 4(a2) and 4(c2) (i.e. an increase in coherence). Additionally, it can be seen that there are fewer pronounced highs and lows, on the whole, in the uncoupled system than in the coupled system; in Fig. 2 highs and lows are sporadically interrupted by stints of somewhat suppressed fluctuations about the fixed point. This in turn contributes to the lower overall amplification of the uncoupled system in comparison to that of the coupled system.

III.3 Transcriptional/translational delays can lead to out-of-phase oscillations despite Delta-Notch coupling

We demonstrated in the previous sections that two cells with slightly disparate oscillatory frequencies could synchronise when coupled via Delta-Notch signalling. As the strength of the Delta-Notch coupling is increased, a common frequency is converged upon and the phase lag between the two cells is reduced. Although this is characteristic of many model parameter sets, it is not always the case. As was noted also in Lewis (2003) and Momiji and Monk (2009) (in the purely deterministic setting), neighbouring cells can be coupled in such a way that they oscillate in anti-phase with one another. This type of behaviour is facilitated by the delays associated with translation and/or transcription.

Fig. 6 demonstrates that for certain parameter sets, as the coupling is increased, oscillations in two cells can converge on a common frequency but the phase lag between the two cells can approach values closer to $\phi=\pi$ than $\phi=0$. That is, the cells tend towards oscillating in anti-phase with one another. Clearly, this is suboptimal if these cellular oscillations are to be used as a segmentation clock.

To understand why this should happen, we monitor the dominant frequency and the associated phase difference between the two cells as the transcriptional/translational delays are varied (see Fig. 7). For the purposes of this analysis, the two cells are taken to be identical. In this case, the two cells are guaranteed to share a peak oscillatory frequency $\omega_{\mathrm{max}}^{(1)}=\omega_{\mathrm{max}}^{(2)}=\omega_{\mathrm{max}}$. We find that whether the cells oscillate in or out-of phase is determined by an interplay between the delays and the dominant frequency of oscillation.

One observes that the phase lag between the two cells switches between [math] and $\pi$ (and vice versa) when the total delta delay time $\tau^{\mathrm{tot}}=\tau^{(m)d}+\tau^{(p)d}$ reaches certain values $\tau^{\mathrm{tot}}_{i}$, where $i=1,2,3,\dots$ (we notice a similar effect when other pairs of delays are varied). In Figure 7, the times at which the switches occur are separated by a regular time interval such that $\tau^{\mathrm{tot}}_{i}-\tau^{\mathrm{tot}}_{i-1}\approx 30\mathrm{mins}$. This interval is roughly equal to half of the typical period of the oscillations, $T/2=\pi/\omega_{\mathrm{max}}$, which varies within around $\sim 10\%$ of the mean value $\sim 30\mathrm{mins}$.

In order to gain some intuition for what this means, we suppose that for a particular set of delays, the two cells oscillate in phase. If we add an additional $\pi/\omega_{\mathrm{max}}$ to the total delay time, the ‘signal’ from one cell to its neighbour is delayed by half a cycle. If the cells would oscillate in phase without this additional delay, it stands to reason that they would oscillate in anti-phase given the additional delay – there would be no difference from the point of view of either cell (assuming that the frequency of oscillation does not vary greatly with the changing delay times). A similar argument holds true for the addition of $2\pi/\omega_{\mathrm{max}}$ to $\tau^{\mathrm{tot}}$ – in this case the phase difference between the two cells ought to be unchanged. One caveat to this reasoning is that the frequency of the oscillations cannot be changed too drastically by the variation in time delay.

We conclude that while the time delays associated with translation and transcription are crucial for the persistence of the cycles which constitute the cell clock, they are somewhat of a double-edged sword. Depending on the interaction between the delays and the internal clocks of each cell, the cells may oscillate in anti-phase with one another. As a result, the precise nature of the transcriptional/translational delays is of great importance with regards to the proper synchronisation of the segmentation clock.

III.4 Transcriptional/translational delays can disrupt global oscillations in chains of cells and give rise to waves of gene expression

In this section, we discuss how the preceding analysis can be extended to a chain of Delta-Notch-coupled cells. We demonstrate that synchronised noisy oscillations in the two-cell system can correspond to global oscillations in a chain. We also explore behaviours other than global oscillations which can occur as a result of delays; namely, the emergence of noise-induced waves.

A distinction between oscillations of the deterministic trajectory and purely noise-induced oscillations was made in Section III.1. In a similar way, one finds that waves can manifest in the individual-based system when they do not in the deterministic system. Such ‘stochastic waves’ or ‘quasi-waves’ have been found previously in theoretical models of individual-based systems with long-range interaction Biancalani et al. (2011). Here however, the stochastic waves arise due to a combination of the non-local dependence of the reaction rates (due to Delta-Notch signalling) and the transcriptional/translational delays. Conversely, waves of gene expression have been studied previously in chains of coupled genetic oscillators Momiji and Monk (2009) but this was done in the context of deterministic equations which ignored intrinsic noise.

We mention the emergence of waves here not so much as an explanation for the travelling waves which are seen in the PSM (these are most likely due to a variation in translatonal/transcriptional delay along the anterior-posterior axis Keskin et al. (2018); Ay et al. (2014)), but as an illustration of the different kinds of undesirable behaviour which can arise when cells oscillate out of phase with one another. As such, the results presented in section are exploratory and not necessarily an attempt to recreate any (as yet) observed phenomena.

We find that for sets of parameters where one would observe oscillations in anti-phase in the two-cell system, one finds waves of gene expression in an extended chain of cells. For parameter sets where the cycles in the two-cell system oscillate in unison, one observes global in-phase oscillations in gene expression.

For an extended chain of cells, we define the power spectrum

[TABLE]

where we have used the discrete Fourier transform with respect to the cell number $j$ defined by $\tilde{f}_{k}=\frac{1}{\sqrt{L}}\sum_{j}e^{ijk}f_{j}$, where $L$ is the number of cells in the chain. Details of the calculation of this power spectrum are given in the Supplement (Section S2 B).

In a chain of coupled cells, global in-phase oscillations are characterised by a peak in the power spectrum $P_{k}^{\alpha}\left(\omega\right)$ at spatial wavenumber $k=0$ and non-zero oscillatory frequency $\omega$. Such a power spectrum is shown in Fig. 8(a), and an example of the corresponding behaviour in a chain of cells is demonstrated in Fig. 8(b), where one observes that the peaks and troughs in one cell tend to align with those in the neighbouring cells. That the cells are indeed in phase with one another is verified by the phase lag (see inset of Fig. 8(a)) which is equal to zero, regardless of cell separation.

Stochastic waves, on the other hand, are characterised by a peak in the Fourier power spectrum at non-zero values of both the spatial wavenumber $k$ and the temporal frequency $\omega$. An example of such a power spectrum is given in Fig. 9(a). To validate the claim that this peak in the power spectrum is indicative of travelling waves, we note that the phase difference between cells varies linearly with cell separation, as is shown in the inset of Fig. 9(a). We stress that the coupling between the cells here is biased in one direction, which breaks the symmetry of the system, allowing waves to travel. Travelling waves of gene expression have been observed in experimental systems other than the PSM Wenden et al. (2012); Ukai et al. (2012).

For symmetric coupling however, one instead observes standing waves of gene expression, where alternate cells oscillate in antiphase, as shown in Fig. 10. In this particular case, the phase lag between any pair of adjacent cells (at the peak frequency) is $\pi$, as can be seen in the inset of Fig. 10(a). This is rather reminiscent of the on-off chequerboard patterns associated with neural differentiation and lateral inhibition Collier et al. (1996).

Examples of travelling and standing stochastic waves in a chain of coupled cells are shown in Figs. 9(b) and 10(b) respectively. There is a clear qualitative distinction between the two. In Fig. 9(b), peaks and troughs in one cell gradually travel in the positive $j$ direction as time goes by, indicating a travelling wave. In Fig. 10(b) however, the peaks in one cell tend to line up with the troughs of the neighbouring cells and visa versa and there is no clear direction of travel.

IV Summary and discussion

The process of somite segmentation poses a complex many-faceted problem for theorists and experimentalists alike and remains an area of active inquiry. Much is left to be discovered about the precise nature of the role of each of the genes involved with the somite segmentation clock and their interactions with external signalling factors in the embryo.

In this paper, we aim to have provided some insight into the role that Delta-Notch signalling plays in not only aligning the frequencies of oscillation of cyclic gene expression in neighbouring cells but also in reducing phase lag, in improving the coherence of oscillations and in increasing their amplitude, so as to produce a robust and reliable segmentation clock. We also explored the role that intrinsic noise plays in the system; counter to intuition, it can actually promote persistent cycles, rather than obscure them. Further, we discussed how the delays involved in the transcription and translation processes can act to promote oscillations but can also result in neighbouring cells oscillating out-of-phase with one another. We examined how this resulted from an interplay between the dominant frequency of oscillation in the cells with the aggregate time-delay. We went on to show how asynchronous behaviour in a two-cell system corresponds to waves of gene expression in a chain of cells.

In a recent work Keskin et al. (2018), the gene expression noise in the PSM of the zebrafish was analysed. This was done by using smFISH microscopy techniques to count the numbers of discrete RNA molecules in individual cells. The statistical discrepancies between the gene expression in sets of cells which were supposedly synchronised was then evaluated (using well-known techniques Elowitz et al. (2002); Swain et al. (2002)) and termed ‘expression noise’. It was found that this expression noise increased when mutations in both the DeltaC and DeltaD genes were introduced, reducing the efficacy of the Delta-Notch coupling. In our work, we have shown that increased Delta-Notch signalling strength increases the degree of synchronisation and the coherence of noisy oscillations, which in turn reduces the discrepancy between the cycles in coupled cells. This would appear to be very much in keeping with the aforementioned experimental findings.

On a more general note, intrinsic noise is often assumed to be a destabilising influence on cycles and to be the source of the discrepancy between the expression in two otherwise-equal cells. But, because of the complexity of the gene regulatory network and the nature of the coupling between cells in the PSM, the intrinsic noise can actually give rise to correlated sustained oscillations in neighbouring cells– a behaviour one might normally associate with an extrinsic influence. This rather blurs the line between the what might be considered the signatures of intrinsic, extrinsic and ‘expression’ noise in experimental data. As a result the utmost care must be taken to identify sources of correlation between cells, other than common external influence, if one is to truly discern the fingerprints of intrinsic noise in the data from those of extrinsic stochasticity.

Author contributions

JWB designed the study, contributed to discussions guiding the work, carried out mathematical calculations, performed simulations and wrote the manuscript. TG designed the study, contributed to discussions guiding the work and wrote the manuscript.

Funding statement

JWB thanks the Engineering and Physical Sciences Research Council (EPSRC) for funding (PhD studentship, EP/N509565/1). TG acknowledges partial financial support from the Maria de Maeztu Program for Units of Excellence in R&D (MDM-2017-0711). We declare no competing interests.

S1 Details of the model gene-regulatory system

In our analysis, we capture the intrinsic noise (which comes about due to the stochastic nature of the transcription/translation processes and finite particle numbers) using an individual-based model. In this model, we suppose that individual protein and mRNA molecules can be created or annihilated with certain probabilities per unit time, which may depend upon the various numbers of proteins/mRNAs currently in existence. There may also be delays between the initialisation of a creation/annihilation event and its completion.

The dynamics, depicted in Fig. 1, are given more precisely by the following set of reactions,

[TABLE]

where $M_{d}$ denotes a molecule of delta mRNA and $P_{h}$ denotes a molecule of Hes/Her protein, etc. The single arrows $\stackrel{{\scriptstyle R}}{{\longrightarrow}}$ indicate that the reaction occurs without delay with a per capita rate $R$. The double arrows $\stackrel{{\scriptstyle[}}{{\tau}}]{R}{\Longrightarrow}$ indicate a delayed reaction with per capita rate $R$ and delay $\tau$. These equations are to be interpreted in the usual way using mass action kinetics Érdi and Tóth (1989). For example, the first reaction is triggered with rate $a_{d}n^{(m)d}$, where $n^{(m)d}$ is the number of $M_{d}$-particles in the system. The effect of such a reaction is realised $\tau^{(p)d}$ units of time later, and results in the addition of a $P_{d}$-particle to the system.

The individual-based dynamics summarised by Eq. (S1) can be approximated by a set of stochastic differential equations (SDEs), which are given in Section S2.3.

The model parameters $a_{\rm h},a_{\rm d}$, $b_{\rm h},b_{\rm d}$ and $c_{\rm h},c_{\rm d}$ are positive rate constants. The composite Hill functions $f^{h}$ and $f^{d}$, encapsulating the activation/inhibition of mRNA production by the various protein concentrations are given by

[TABLE]

where $r_{0}^{\alpha},r_{\rm d}^{\alpha},r_{\rm h}^{\alpha},r_{\rm hd}^{\alpha}$ and $k_{\alpha}$ are rate constants, $N$ is the system size, $\phi^{\alpha}_{j,t}=\frac{n^{\left(p\right)\alpha}_{j,t}}{n^{\left(p\right)\alpha}_{0}}$, and where $d^{(\pm)}$ are positive constants such that $\frac{1}{2}\left(d^{(+)}+d^{(-)}\right)=1$; $n^{\left(p\right)\alpha}_{0}$ are reference values for the protein levels. Superscript or subscript indices $\alpha$ are placeholders, representing the cases $\alpha\in\{{\rm h,d}\}$. This reaction scheme is similar to the one used in Lewis (2003).

The different terms in the function $f^{\alpha}_{j}\left(\mathbf{n}_{t}\right)$ correspond to autonomous activation ($r^{\alpha}_{0}$), activation by delta only ($r^{\alpha}_{\mathrm{d}}$), inhibition by hes/her only ($r^{\alpha}_{\mathrm{h}}$) and mixed response to hes/her and delta $r^{\alpha}_{\mathrm{hd}}$. We constrain the parameters $r_{0}^{\alpha},r_{\rm d}^{\alpha},r_{\rm h}^{\alpha},r_{\rm hd}^{\alpha}$ always to sum to unity i.e. $r^{\alpha}_{0}+r^{\alpha}_{\mathrm{d}}+r^{\alpha}_{\mathrm{h}}+r^{\alpha}_{\mathrm{hd}}=1$. The expression inside the square brackets in Eq. (S2) can therefore range between [math] and $1$ dynamically, depending on the concentrations of mRNA or protein and on the parameters $r_{0}^{\alpha},r_{\rm d}^{\alpha},r_{\rm h}^{\alpha},r_{\rm hd}^{\alpha}$. As a result, the typical mRNA birth rate is characterised by $k_{\alpha}N$.

In the main text, we examine the change in the behaviour of coupled genetic oscillators as we vary the ‘coupling strength’. In order to isolate the effect of the Delta-Notch signalling from the hes/her auto-repression, we always set $r^{\mathrm{h}}_{0}=r^{\mathrm{h}}_{\mathrm{d}}=0$ and vary $r^{\mathrm{h}}_{\mathrm{hd}}$ and $r^{\mathrm{h}}_{\mathrm{h}}$ such that $r^{\mathrm{h}}_{\mathrm{h}}+r^{\mathrm{h}}_{\mathrm{hd}}=1$. As $r^{\mathrm{h}}_{\mathrm{hd}}$ increases, so does the coupling strength but the role of hes/her remains the same. For zero coupling $r^{\mathrm{h}}_{\mathrm{hd}}=0$. For maximal coupling $r^{\mathrm{h}}_{\mathrm{hd}}=1$. We keep the typical mRNA production rates $k_{\alpha}$ constant.

In Sections III A, III B 1, III B 2 and III C, since we only consider a 2-cell model,     $\frac{1}{2}\left(d^{(-)}\phi^{\mathrm{d}}_{j-1,t}+d^{(+)}\phi^{\mathrm{d}}_{j+1,t}\right)=\phi^{\mathrm{d}}_{j^{\prime},t}$ where $j^{\prime}=2$ if $j=1$ and vice versa.

S2 Quantification of stochastic fluctuations in systems with delays and non-local reaction rates

In this section, we derive the analytical results for the quantification of the stochastic fluctuations about the fixed point of a delay system with non-local reaction rates (Delta-Notch signalling). First, using a path integral approach, we derive expressions for the effective stochastic differential equations (SDEs) Gardiner (2009) which approximate the individual-based dynamics in the limit of large system-size $N$. We then use the linear-noise approximation to obtain an expression for the correlators of fluctuations about the deterministic fixed-point of the system. This procedure is similar to that used in Brett and Galla (2013, 2014); Baron and Galla (2019). These correlators are the basis for the theory results presented in the main text. Finally, we detail how this analysis can be used in conjunction with the model detailed in Section S1 to produce the results in the main text.

S2.1 Path-integral approach

We begin our analysis by defining a stochastic process in terms of the scaled variables $q^{\alpha}_{j,t}=n^{\alpha}_{j,t}/N$ (we use a more compact notation for the numbers of molecules $n^{\alpha}_{j,t}\equiv n^{\alpha}_{j}(t)$ than in the main text). Here, $N$ (defined in Section S1) characterises the typical number of particles per cell and is sometimes referred to as the system size Van Kampen (1992). To simplify matters, we discretise time in steps of length $\Delta$. The continuum limit is later recovered by taking $\Delta\to 0$. When reactions occur at a site $j$, particles may be created or annihilated immediately and/or at one other future time $t+\tau$. We define the number of reactions of type $r$ which occur in the time interval $t$ to $t+\Delta$ and have an associated delayed effect at $t+\tau$ by $k_{j,r,\tau,t}$. The number of particles of type $\alpha$ [where here $\alpha\in\{(p)d,(p)h,(m)d,(m)h\}$] which are immediately created/annihilated in such a reaction is denoted $\nu_{r}^{\alpha}$ and the number which are created/annihilated at the later time $t+\tau$ is denoted $w_{\tau,r}^{\alpha}$.

The stochastic process can therefore be written as

[TABLE]

We wish to approximate this process with a set of stochastic differential equations. An elegant way to obtain this approximation is to start from a path-integral representation for the process, and then to perform an expansion in inverse powers of the system size (similar to that used by van Kampen Van Kampen (1992)) within this representation. Using a path-integral approach, as opposed to a master equation, avoids the complications which arise due to the non-Markovian nature of the dynamics.

We can write an expression for the probability of observing the system in a configuration $\mathbf{n}$ (which represents set of particle numbers of all types and locations) at time $t$ as follows

[TABLE]

where $\mathcal{P}\left(\{q^{\alpha}_{j,t^{\prime}}\}\right)$ is the probability density of observing a particular trajectory (realisation of the system) $\{q^{\alpha}_{j,t^{\prime}}\}$ and $\delta(\cdot)$ is the Dirac delta-function. Eq. (S4) is merely a statement that the probability of the system being in state $\mathbf{n}$ at time $t$ is equal to the sum of the probabilities of observing any one of a set of paths which lead to the system being in state $\mathbf{n}$ at time $t$.

Constraining the trajectory to obey Eq. (S3) and then rewriting the Dirac delta functions as integrals of complex exponentials, one obtains

[TABLE]

where $P\left(\{k\}\right)$ is the probability of observing a particular set of creation/annihilation events $\{k\}$.

We presume that each reaction event is independent such that

[TABLE]

Once we specify the exact probability distributions for the variables $\{k\}$, we can evaluate the sums over $\{k\}$ in Eq. (S5). Let $W_{j,r,\tau,t}$ be the probability per unit time squared (or ‘rate’) of a reaction of type $r$ occurring at position $j$ at a time between $t$ and $t+\Delta$ with associated delay between $\tau$ and $\tau+\Delta$. These rates may be dependent on the numbers of particles $\mathbf{n}$. We presume that the numbers of reactions triggered in the interval $t$ to $t+\Delta$, $k_{j,r,\tau,t}$, are Poisson random variables with mean $W_{j,r,\tau,t}\Delta^{2}$. This involves the approximation that the reaction rates do not change over the course of the small time step $\Delta$. This is similar to the approximation made for the so-called tau-leaping variant of the Gillespie algorithm Gillespie (2001). Making these assumptions, the distributions of the $\{k\}$ are given by

[TABLE]

We note that in our system, the local reaction rate $W_{j,r,\tau,t}$ can depend upon the number of particles in the adjacent cells. The sums over $\{k\}$ in Eq. (S5) can be evaluated by observing that

[TABLE]

We finally arrive at the following expression for the probability distribution $P\left(\mathbf{n},t\right)$

[TABLE]

Expanding the exponentials in Eq. (S9) in powers of $1/N$ and truncating the series at next-to-leading order, one obtains the following expression,

[TABLE]

where $\delta_{l,l^{\prime}}$ is the Kronecker delta. Eq. (S10) is similar to the Martin-Siggia-Rose-Janssen-de Dominicis (MSRJD) functional integral Altland and Simons (2010). This result can be compared with the path integral for an SDE of an appropriate form. The corresponding path integral expression for the general SDE

[TABLE]

is given by Brett and Galla (2013, 2014); Baron and Galla (2019)

[TABLE]

From Eq. (S10) one can then read off the following effective SDEs, which are good approximations of the dynamics for large but finite $N$, by taking the limit $\Delta\to 0$

[TABLE]

where $W_{j,r,t}=\int_{0}^{\infty}W_{j,r,\tau,t}d\tau$, and where the correlators of the stochastic noise variables are given by

[TABLE]

S2.2 Linear-noise approximation and calculation of the power spectra and phase lags

We presume that the reaction rates can be decomposed as follows $W_{j,r,\tau,t}=W_{j,r,t}K_{r}\left(\tau\right)$. That is, the delay time $\tau$ is drawn from a distribution $K_{r}\left(\tau\right)$ which is independent of $j$, $t$ and $n^{\alpha}_{j,t}$. Eq. (S13) then becomes

[TABLE]

If we consider small fluctuations about the fixed point of the deterministic system $\delta^{\alpha}_{j,t}=n^{\alpha}_{j,t}-\bar{n}^{\alpha}$ [as in Eq. (1)], Eq. (S15) may be approximated by

[TABLE]

where

[TABLE]

and

[TABLE]

where $\theta\left(\tau\right)$ is the Heaviside function. We note that for the systems studied in the main text, $J^{\alpha,\alpha^{\prime}}_{j,j^{\prime}}$ and $L^{\alpha,\alpha^{\prime}}_{j,j^{\prime},\tau}$ are non-zero for $j\neq j^{\prime}$ due to the coupling of adjacent cells through Delta-Notch signalling.

Crucially, as a part of this approximation, we now neglect fluctuations about the fixed point of the system in the evaluation of the correlators Eq. (S14). That is, we evaluate $\langle\xi^{\alpha}_{j,t}\xi^{\alpha^{\prime}}_{j^{\prime},t^{\prime}}\rangle$ at $n^{\alpha}_{j,t}=\bar{n}^{\alpha}$. As result, what was multiplicative noise is now treated as additive, simplifying the calculation.

Carrying out a temporal Fourier transform in Eq. (S16), one obtains

[TABLE]

Eq. (S19) can be rewritten in matrix form as

[TABLE]

where the different elements of the vector correspond to the different types of particle at the different cell sites. Writing the correlation matrix of the noise variables as $\underline{\underline{\Sigma}}_{\omega}=\langle\hat{{\underline{\xi}}}_{\omega}\hat{{\underline{\xi}}}^{\dagger}_{\omega}\rangle$ one finally obtains the following result for the matrix of the correlators of the fluctuations

[TABLE]

The diagonal elements of this matrix correspond to the power spectrum of fluctuations in Eq. (2), i.e. $P^{\alpha}_{j}\left(\omega\right)=C^{\alpha,\alpha}_{j,j}\left(\omega\right)$. The off-diagonal elements allow one to calculate the phase lag as stated in Eq. (3), that is $C^{\alpha,\alpha^{\prime}}_{j,j^{\prime}}\left(\omega\right)=\langle\hat{\delta}^{\alpha}_{j,\omega}\hat{\delta}^{\alpha^{\prime}\star}_{j^{\prime},\omega}\rangle$.

When the problem is translationally invariant, as is the case when the system parameters are the same in all cells, $J^{\alpha,\alpha^{\prime}}_{j,j^{\prime}}$ and $L^{\alpha,\alpha^{\prime}}_{j,j^{\prime},\tau}$ are functions of $j-j^{\prime}$ only. One can then further simplify matters by carrying out a Fourier transform with respect to position as well. One obtains

[TABLE]

and from this

[TABLE]

where the different elements of the vector now correspond only to the different species. We then arrive at a similar expression to Eq. (S21), but where the matrix dimension is reduced by a factor of $L$ (the number of cells)

[TABLE]

Here, the diagonal elements correspond to the power spectra in Eq. (6) and the off-diagonal elements allow one to calculate the phase lag using Eq. (3) of the main text.

S2.3 Application to the model gene regulatory network

Using Eqs. (S13) and (S14), the individual-based system given by Eq. (S1) can be approximated by stochastic differential equations of the form

[TABLE]

The quantities involved in these equations are given in Section S1, with the exception of the Gaussian noise terms $\xi^{\alpha}_{j,t}$, which encapsulate the intrinsic noise due to the stochastic and individual-based nature of the gene regulatory system. The deterministic trajectories $\bar{n}^{\alpha}_{j,t}$, examples of which are given in Figs. 2 and 3, are found by setting the noise terms $\xi^{\alpha}_{j,t}$ in Eq. (S25) to zero and numerically integrating the resulting ordinary differential equations. Setting the noise terms to zero effectively approximates the system as being infinitely large, so that any fluctuations are negligible in comparison to the mean particle numbers. This is not appropriate in cases in which intrinsic noise significantly affects the dynamics (as is the case for the examples we look at). This is exemplified by the stark disagreement between the deterministic and individual-based simulations in Figs. 2 and 3.

When the deterministic (infinite) system has reached a fixed point, the noise terms $\xi^{\alpha}_{j,t}$ in Eq. (S25) can be taken to have the following correlators within the linear-noise approximation,

[TABLE]

where barred quantities are evaluated at the deterministic fixed point (i.e. the solution to Eqs. (S25) with the noise terms and the time derivatives set to zero). All inter-species cross-correlators are zero, due to the fact that no one reaction gives rise to the production/annihilation of two different species of particle (see Eq. S1).

In order to evaluate the power spectrum of fluctuations or to find the phase lag, one first finds the matrix $\underline{\underline{M}}_{\omega}$ or $\underline{\underline{M}}_{k,\omega}$, defined in Eqs. (S20) or (S23) respectively, by performing the linear-noise approximation on Eq. (S25), as detailed in Section S2.2. Using the correlators in Eq. (S26), one is then able to evaluate the matrix of correlators $\underline{\underline{C}}\left(\omega\right)$ or $\underline{\underline{C}}\left(k,\omega\right)$, defined in Eqs. (S21) or (S24) respectively, which contains all the information one needs to find the power spectra of fluctuations and/or the phase lag (as is also described in Section S2.2).

As an example, we provide expressions that allow one to calculate the matrix $\underline{\underline{M}}_{k,\omega}$ in the case that the system is spatially homogeneous. In Eq. (S22), the matrix $\tilde{J}_{k}^{\alpha\alpha^{\prime}}$ is given by

[TABLE]

and the matrix $\hat{\tilde{L}}^{\alpha\alpha^{\prime}}_{k,\omega}$ is given by

[TABLE]

where the entries of the matrix are in the order $\alpha=\{\mathrm{(p)h},\,\,\mathrm{(p)d},\,\,\mathrm{(m)h},\,\,\mathrm{(m)d}\}$ and where we have

[TABLE]

S3 Linear stability analysis

We have derived expressions for the power spectrum of fluctuations in the individual-based system about the deterministic fixed point. Such an analysis presumes the stability of this fixed point. To determine whether or not such a stable fixed point exists for a particular parameter set, we perform a linear stability analysis of the deterministic system.

We begin with the linearised expression for the time-evolution of small deviations about the deterministic fixed point Eq. (S16), but with the noise term removed

[TABLE]

For simplicity, we now assume that the delay kernels are all Dirac-delta functions [i.e. $K_{r}\left(\tau\right)=\delta\left(\tau-\tau_{r}\right)$] such that

[TABLE]

where

[TABLE]

We have further supposed that the number of particles produced by a delayed reaction $w^{\alpha}_{r}$ is independent of the delay time $\tau_{r}$.

We then suppose that the deviation away from the fixed point evolves in an exponential way. That is we propose the ansatz

[TABLE]

Upon substitution into Eq. (S30), one obtains

[TABLE]

We note that if we had not assumed a delta function for the delay kernel $K_{r}\left(\tau\right)$, the exponential factors in Eq. (S35) would instead be replaced by $\mathcal{L}_{\tau}\left\{K_{r}\left(\tau\right)\right\}\left(-\lambda\right)$, the Laplace transform of the delay kernel evaluated at $-\lambda$. Eq. (S34) can be rewritten more succinctly as the following matrix equation

[TABLE]

In order for the solutions to this equation to be non-trivial, the determinant of the object multiplying $\underline{\delta}^{(0)}$ in Eq. (S35) must be equal to zero. This gives rise to an ‘eigenvalue’ equation for $\lambda$. Due to the exponential terms involving $\lambda$ which arise from the delays, this equation is not analytically tractable and must be solved numerically.

Typically, one finds many possible solutions $\lambda$ to the eigenvalue equation for any one set of system parameters, and thus the full solution to Eq. (S30) is a linear combination of solutions of the form in Eq. (S33). If any one of the eigenvalues $\lambda$ has a positive real part the fixed point is unstable. Else, we say that it is (linearly) stable. In order to produce the dotted lines in Figs. 5 and 6, one finds the eigenvalues $\lambda$ as a function of the system parameters. The instability line divides regions in parameter space for which all the eigenvalues have negative real parts from those regions for which some eigenvalues have positive real parts. The imaginary part of the least stable eigenvalue for each parameter set is plotted as a green line in Figs. 5a and 6a.

Once one has found the eigenvalues $\lambda$, one can then also go on to solve (numerically) for the eigenvectors $\underline{\delta}^{(0)}$. These yield information about the relative amplitudes and phases of the various components. The complex phase difference between elements of the eigenvector that correspond to the same species in opposite cells is shown as a green line in Figs. 5b and 6b. The inter-cell phase difference turns out to be the same for all particle types.