Intrinsic noise, Delta-Notch signalling and delayed reactions promote sustained, coherent, synchronised oscillations in the presomitic mesoderm
IIntrinsic noise, Delta-Notch signalling and delayed reactions promote sustained,coherent, synchronised oscillations in the presomitic mesoderm
Joseph W. Baron ∗ and Tobias Galla
1, 2, † Theoretical Physics, School of Physics and Astronomy,The University of Manchester, Manchester M13 9PL, United Kingdom Instituto de F´ısica Interdisciplinar y Sistemas Complejos IFISC (CSIC-UIB), 07122 Palma de Mallorca, Spain
Using a stochastic individual-based modelling approach, we examine the role that Delta-Notchsignalling plays in the regulation of a robust and reliable somite segmentation clock. We find that notonly can Delta-Notch signalling synchronise noisy cycles of gene expression in adjacent cells in thepresomitic 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 anddemonstrate how intrinsic noise can play an active role in promoting sustained oscillations, givingrise to noise-induced quasi-cycles. Finally, we explore how translational/transcriptional delays canresult in the cycles in neighbouring cells oscillating in anti-phase and we study how this effect relatesto the propagation of noise-induced stochastic waves.
I. INTRODUCTION
In developing vertebrates and cephalochordates, as theembryo forms and extends pairs of blocks of mesodermalprogenitor cells assemble, bilaterally flanking the noto-chord [1]. These blocks, termed somites, eventually go onto form vertebrae and ribs after further cellular differenti-ation. The somites are constructed pair-by-pair, anteriorto posterior, in a rythmic and sequential manner as thetailbud extends away from the rostral end of the embryo.They are formed from cells originating in the presomiticmesoderm (PSM). Such cells are produced continually bythe tailbud as the abdomen elongates [2, 3].The process of somite segmentation has been of inter-est to experimentalists and theorists working in the fieldof developmental biology for some decades; it provides afascinating case study where one can directly examine thelink between microscopic gene regulatory systems oper-ating in individual cells and macroscopic developmentalprocesses. The prevailing theoretical framework for un-derstanding the process was put forward by Cooke andZeeman in 1976 [4] and is termed the ‘clock-wavefront’model. This model proposes that the cells in the pre-somitic mesoderm each possess an internal cyclic ‘clock’which is synchronised between the cells. Additionally, awavefront propagates through the PSM as the embryogrows. As the wavefront encounters cells, it interactswith them differently depending on the current state ofeach internal cellular clock. This interaction causes thecells to change their adhesive and migratory properties.The temporal periodicity of the cell cycles is thus con-verted into the spatial periodicity of the somites.Considerable experimental and theoretical effort hasbeen expended in order to identify the genetic oscillatorsthat constitute the putative somite segmentation ‘clock’and a good amount of progress has been made. In certain ∗ [email protected] † [email protected] model organisms, such as the mouse and the zebrafish,so-called ‘knockdown’/‘knockout’ experiments have iden-tified genes which when mutated give rise to defects inthe formation of somites and, consequently, the verte-brae [5–11]. Gradients of FGF (fibroblast growth factor)or Wnt protein, which are produced in the tailbud, arethought to constitute the moving wavefront; transientloss or increase in these substances can alter the localsomite length [12, 13]. Genes such as hes in the mouse[14] and her in the zebrafish [15] are thought to be theprimary cyclic genes which act as clocks. These genes areboth targets of the Notch signalling pathway. It has alsobeen shown in experiments that Delta-Notch signalling isa vital component in synchronising oscillations [14, 16–18].In order for the oscillations in the expression of the hes/her genes to constitute a viable segmentation clockfor the clock-wavefront model, the oscillations must sat-isfy several criteria: (1) The oscillations in gene expres-sion 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 bea clear dominant frequency). (4) The oscillations musthave sizeable enough an amplitude so as not be indistin-guishable from background ‘noise’. Mathematical modelsof the gene regulatory system have shown that Delta-Notch signalling can indeed synchronise (align the fre-quencies) of oscillations in neighbouring cells with intrin-sically differing cellular clocks [19]. That is, it has beenshown that Delta-Notch signalling is responsible for sat-isfying condition (1), but relatively little discussion hasbeen dedicated to the latter 3 conditions (phase, coher-ence and amplitude).So that one might analyse the degree to which thecriteria above are satisfied, one must take into accountstochastic (random) effects in the system, especially withregards to point (3). The gene regulatory systems inquestion are inherently noisy in nature [20–26]. Thisis due in part to the stochastic nature of the produc-tion/decay events of individual proteins and/or mRNA a r X i v : . [ q - b i o . S C ] N ov molecules and the fact that there are finite numbers ofthese molecules in any one cell. Noise of this origin istermed intrinsic in the literature [27, 28]. On the otherhand, gene expression is also influenced by the concentra-tions, locations and states of molecules such as regulatoryproteins or polymerases which can affect the global activ-ity in a single cell but can vary between cells. Noise aris-ing from fluctuations in the properties of such moleculesis referred to as extrinsic [27, 28].The role of noise has largely been disregarded in pre-vious theoretical work on the somite segmentation clock[29–31], or has often been treated only as an externalinfluence rather than as an aspect intrinsic to the trans-lation and transcription processes [32–34] (an exceptioncan be found in [35], where simulations involving intrin-sic noise were performed). For example, some works haveconsidered the binding and unbinding of repressor proteinto the DNA binding site as a stochastic process [19, 36] sothat the rates of transcription are themselves stochasticvariables. This source of noise is taken into account inthe context of deterministic evolution equations - the in-trinsic stochasticity of the transcription and translationevents is not accounted for. In this paper however, westudy the effect of intrinsic stochasticity by treating theproduction and the decay of individual molecules as ran-dom processes, following [20, 37–41]. These events maybe subject to delays arising from the finite time taken forthe translation/transcription processes.Using an individual-based mathematical model captur-ing the intrinsic noise in the system, we demonstrate that,perhaps counter to intuition, intrinsic stochasticity canbe a proactive force in promoting cellular oscillations.We study how these noisy oscillations in neighbouringcells are affected by different levels of the strength ofDelta-Notch signalling. We are able to show that, undercertain conditions, Delta-Notch signalling not only actsto align the frequencies of oscillations in neighbouringcells; it can also reduce the phase lag, reduce the rangeof dominant oscillatory frequencies (i.e., it can make theoscillations more coherent) and it can increase the am-plitude of oscillations (also noted in [17, 35]). The com-bination of these effects indicates that Delta-Notch sig-nalling can contribute to satisfying points (1)-(4) above.We also discuss circumstances under which pairs of cellsmay oscillate out-of-phase, despite Delta-Notch coupling.We explore how this is related to waves and/or oscillat-ing chequerboard patterns of gene expression in extendedchains of cells. II. METHODSA. Model definition
Oscillations in the expression of genes (or ‘pulsing dy-namics’) is a well-documented phenomenon responsiblefor many cell functions and broader biological processes[42, 43]. Such oscillations in the expression of certain genes are thought to constitute the biological ‘clock’ inthe clock-wavefront model [4] of somite segmentation[2, 3]. The genes in question are known to be affectedby Notch signalling. For the purposes of our theoreticaltreatment it is not necessary to consider the full complex-ity of the Notch signalling pathway [44, 45] or even thefull network of interacting genes involved with the somitesegmentation process [46]; one can use a simplified modelto highlight the salient features and analyse their causes.So-called ‘knockdown’/‘knockout’ experiments [3] sug-gest that the most relevant genes for the regulation ofthe ‘clock’ are delta (or its homologues) and hes in mice,and her in zebrafish [2]. It has been shown previously [31]that a two-gene model involving only hes/her and delta is sufficient for the emergence of cycles. We therefore alsoadopt a reduced two-gene model, as this will be sufficientto highlight the effects with which we are concerned. Thereduced system is depicted schematically in Fig. 1 and isdiscussed in more detail in the Supplement (Section S1).For now, we consider a system of two coupled cells as asimple example. We generalise the approach to systemsof greater numbers of cells in Section III D.
FIG. 1. Schematic of the reduced 2-cell gene regulatory sys-tem [19, 31]. Genes ( hes/her and delta ) are transcribedto produce mRNA. In turn, protein is translated from themRNA, which goes on to activate/inhibit further mRNA tran-scription. Both the transcription and translation processestake an amount of time to complete, giving rise to delays τ ( m ) x and τ ( p ) x respectively, where x ∈ { h, d } . Hes/Her pro-tein inhibits local delta transcription. Delta protein acts as aligand to the Notch receptor on the adjacent cell. Notch inturn activates the production of hes/her mRNA. IndividualmRNA and protein molecules degrade (and become inactive)at constant probabilities per unit time. Using n ( t ) to denote the set of all protein and mRNAnumbers in all cells at time t , the reduced gene reg-ulatory system that we consider can be summarisedas follows: hes/her mRNA is transcribed at a rate f h (cid:2) n (cid:0) t − τ ( m ) h (cid:1)(cid:3) . The rate f h (cid:2) n (cid:0) t − τ ( m ) h (cid:1)(cid:3) takes theform of a sum of Hill functions, which reflects the factthat hes/her mRNA production is inhibited by localHes/Her protein and activated by Delta protein in ad-jacent cells. The precise form of this non-linear func-tion is given in the Supplement (Section S1). Every hes/her mRNA molecule is translated into protein at aconstant rate a h . Because the transcription and transla-tion processes take an amount of time τ ( m ) h and τ ( p ) h to complete respectively, the present rate of produc-tion of mRNA/protein is dependent upon protein/mRNAconcentrations 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. Ina similar way, delta mRNA is transcribed at a rate f d (cid:2) n (cid:0) t − τ ( m ) d (cid:1)(cid:3) and decays at a constant per capitarate c d . Delta protein molecules are produced and decayat the per capita rates a d and b d respectively. Produc-tion of delta mRNA and protein are delayed by the times τ ( m ) d and τ ( p ) d to respectively. Finally, production of delta mRNA is inhibited by local Hes/Her protein con-centration.There are three vital aspects to the processes in thissetup: (1) The model is individual-based – it does nottreat protein/mRNA concentrations as continuous quan-tities (an approximation only valid when populationnumbers are large). The production and degradationof proteins and mRNA are inherently stochastic (ran-dom) processes due to the finite numbers of proteins andmRNA [23, 24] in each cell; this gives rise to noisy dynam-ics [40]. (2) There is a time-delay between the activationof the production of one unit of mRNA/protein and thecompletion of the production process. As a result, therates of production of mRNA/protein at a given timeare dependent on the state of the system in the past. Inthe language of stochastic processes, the dynamics arenon-Markovian (they have memory) [47]. It has been es-tablished that time-delays such as these can encouragethe emergence of temporal oscillations [19, 31, 48]. (3)Due to Delta-Notch signalling, the rate of production of hes/her mRNA in one cell is dependent on the concen-tration of Delta protein in the neighbouring cells. In thissense, there is a non-locality to the reaction rates.The combination of these three aspects of the dynam-ics leads to a unique challenge with respect to theoreticalmodelling. However, we demonstrate in the Supplementthat one can approximate the full individual-based dy-namics of the system with a set of stochastic differen-tial equations (SDEs); these are given in the Supplement(Section S2 C). They take a similar form to the deter-ministic (noiseless) equations given in [19] but includeadditional Gaussian noise terms which take into accountthe intrinsic stochasticity of the system. We emphasizethat the properties of this noise are calculated so as toagree with individual-based simulations of the system;the noise is not added in an ad hoc fashion. The tools weuse to quantify the phenomena induced in the system bynoise are discussed in the next section. B. Analysis of stochastic behaviour
In this work, we will be concerned primarily with thetheoretical analysis of noise-induced cycles of gene ex-pression. These are oscillations which occur in the fullstochastic individual-based model but which are missingin the noiseless deterministic system.The power spectrum of fluctuations about the deter-ministic trajectory will be the main quantitative toolthat we use to analyse the noise-induced phenomena inthe gene regulatory model described in Section II A (andelaborated upon in the Supplement Section S1). We de-note the number of particles of type α in cell j at time t by n αj ( t ). The type of particle indicated by the index α may be mRNA molecules or proteins. The dynamicsof the quantities n αj ( t ) are approximated by the systemof stochastic differential equations (SDEs) in Eq. (S25)(in the Supplement). Further, we write ¯ n αj ( t ) for thenumbers of particles predicted by the corresponding de-terministic model [Eq. (S25), with the noise terms ξ αj ( t )set to zero]. Thus, we define the fluctuations about thedeterministic trajectory as δ αj ( t ) = n αj ( t ) − ¯ n αj ( t ) . (1)The power spectrum of fluctuations is then defined viathe temporal Fourier transform as P αj ( ω ) = (cid:104)| ˆ δ αj ( ω ) | (cid:105) , (2)where the Fourier transform is given by ˆ g ( ω ) = √ π (cid:82) ∞−∞ e iωt g ( t ) dt ; the angular brackets denote the en-semble average over the set of all possible stochastic timecourses of the system. Roughly speaking, the powerspectrum of fluctuations decomposes a time series intoits composite frequencies and quantifies the statisticalcontribution of a particular frequency to the series. Alarge, narrow, unique peak in the power spectrum indi-cates that the frequency at which the peak is located isthe dominant frequency of the time series; a peak of thistype centred on a non-zero frequency therefore charac-terises periodicity. If the deterministic trajectory ¯ n αj ( t )is non-oscillatory, then such a peak in the spectrum ofthe stochastic model indicates noise-induced oscillations.We note that in our analysis the power spectrum is al-ways evaluated in the steady state, i.e. when all transienteffects have decayed sufficiently so as to be negligible.Using a mathematical approach (see Supplement Sec-tions S1, S2 and S3), we are able to predict the powerspectrum of the fluctuations when the deterministic tra-jectory has reached a fixed point (i.e., when ¯ n αj ( t ) = ¯ n αj is constant at long times t ). This theoretical approachgives us a way to identify, without performing time-consuming simulations, what the dominant frequency ofnoise-induced oscillations is and to what extent the otherfrequencies contribute. This analysis is only valid fornoise-induced cycles, i.e. when the oscillations of the de-terministic equations are transient. For the parameterregimes where there are persistent cycles in the deter-ministic equations, we perform a linear stability analysisin order to obtain quantities such as periods of oscillationand phase lags (see Supplement Section S3).It is also possible to quantify the phase lag betweentwo sets of noise-induced cycles in coupled cells using ourtheoretical approach. Following [49, 50], we define thephase lag φ α,α (cid:48) j,j (cid:48) ( ω ) associated with a particular frequency ω between species α in cell j and species α (cid:48) in cell j (cid:48) astan (cid:16) φ α,α (cid:48) j,j (cid:48) ( ω ) (cid:17) = Im (cid:16) (cid:104) ˆ δ αj ( ω ) ˆ δ α (cid:48) (cid:63)j (cid:48) ( ω ) (cid:105) (cid:17) Re (cid:16) (cid:104) ˆ δ αj ( ω ) ˆ δ α (cid:48) (cid:63)j (cid:48) ( ω ) (cid:105) (cid:17) . (3)Phases differing by integer multiples of 2 π are degeneratetherefore, in this paper, we define the phase lag to bein the range [ − π, π ). Notably, the phase lag as definedin Eq. (3) is dependent on ω . As mentioned above, atime series can be thought of as being comprised of asum of cycles with different frequencies ω . The quantity φ α,α (cid:48) j,j (cid:48) ( ω ) is the phase lag between the constituent cyclesof frequency ω in cells j and j (cid:48) . In our analysis, wemay refer to the phase lag of a cell j (with respect toanother cell), which we define as the phase lag at the peakfrequency of the power spectrum of the cell in question, ω ( j )max .Furthermore, following [51], we define the total am-plification of fluctuations for particles of type α in cell j A αj = (cid:90) ∞ P αj ( ω ) dω. (4)This quantity is proportional to the time-averagedsquared displacement of the dynamics from the fixedpoint, i.e., to the variance of the stochastic time series.Finally, again following [51], we also define the coher-ence as the proportion of the power spectrum within afixed range ∆ ω of the peak C αj = 1 A αj (cid:90) ω ( j )max + ∆ ω ω ( j )max − ∆ ω P αj ( ω ) dω. (5)The coherence C αj quantifies how sharply peaked thepower spectrum is - i.e., how narrow the band of dom-inant frequencies is. It has a maximum value of 1 anda minimum value of 0. The choice of ∆ ω is largely im-material provided ∆ ω is small compared to the peak fre-quency. III. RESULTSA. Individual-based models capture noise-driveneffects which are missed by deterministic models
In the systems we are considering, individual cells con-tain of the order of 10-100 mRNA molecules and around 1000 proteins of any one type [23, 24] (see also [19, 52]).As such, the dynamics are inherently noisy. This type ofnoisy dynamics has been observed in experiments mon-itoring gene expression [21–26]. The expression of thesegenes cannot be fully described by the regular, smoothoscillation obtained from integrating deterministic setsof ODEs (as can be seen from Figs. 2 and 3). Instead,a stochastic individual-based model is better suited toqualitatively reproduce the results of experiment.In previous theoretical studies of the somite segmenta-tion clock, noise has mostly been treated as external tothe dynamics [19, 33, 34] or has not been considered at all[31, 53]. An exception to this is [35], in which individual-based simulations were carried out (the consequences ofthe inclusion of intrinsic noise that we discuss here werenot the focus of [35] however). The noise that we usein the present work is rigorously derived as an intrin-sic quality of the stochastic, individual-based dynamicsthemselves. As such, the results of our analysis agreewith fully individual-based simulations of the system (asis demonstrated in Fig. 4).It has previously been observed [19] that noisy exter-nal driving can give rise to sustained oscillations in thesomite segmentation clock. We show that similar noise-induced oscillations can also be produced by the stochas-tic nature of the transcription and translation processesin the segmentation clock themselves. That is, the inclu-sion of intrinsic noise in the theoretical modelling givesrise to sustained noise-induced oscillations which a purelydeterministic 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 [19]), one may observe the noiseless model tend to-wards a stationary fixed–point, only exhibiting transientoscillations which eventually decay. In the correspond-ing individual-based model however, the noise repeat-edly ‘kicks’ the system away from the fixed point. Asa consequence, the oscillations which were transient inthe deterministic model are sustained by the noise. Onethus observes persistent noisy oscillations for the sameparameter set.The oscillations in gene expression observed experi-mentally may very well be noise-induced cycles of thistype. The emergence of such ‘quasi-cycles’ is a well-documented phenomenon which has been previouslystudied in the context of gene-regulatory models [37, 38]as well as in ecological systems [54, 55] and in epidemics[50, 51, 56]. That these are indeed cycles with a periodicnature and not just random white noise is demonstratedby the power spectra of fluctuations (see Fig. 4)– thismatter is discussed further in Sections II B and III B 1.Our theoretical approach to analysing noise-inducedcycles (which is similar to that found in [57–59] and de-tailed in the Supplement) allows us to study the amplifi-cation, synchronisation and coherence of these cycles, asdiscussed in the following sections.
FIG. 2. Oscillations in Hes/Her protein numbers for the sim-plified two-cell gene regulatory system (depicted in Fig. 1)with no Delta-Notch coupling. Panel (a) shows the results ofindividual-based simulations where noisy cycles persist. Thesimulations of the non-linear stochastic model are performedusing a modified version [60] of the Gillespie algorithm [61],which takes into account delayed reactions. This is in contrastto the deterministic trajectory in panel (b), where the oscil-lations are transient and eventually decay to a fixed point. Asmall difference between the delays in either cell gives rise todiffering frequencies of oscillation - the noisy cycles are notsynchronised. This is further illustrated by the power spectraof the cycles in panel (a), which are shown in Fig. 4(a2). Re-ferring to the model specified in Supplement Section S1, therate parameters used here are a α = 4 . b α = 0 . c α = 0 . k α = 3 . α , the system size is N = 10, the referenceprotein levels are n ( p ) h = 4 N and n ( p ) d = 100 N , the delaytimes are τ (p)h = 2, τ (p)d = 5 and τ (m)d = 50 in both cells but τ (m)h = 18 in cell 1 and τ (m)h = 22 in cell 2. The rates asso-ciated with the Hill functions are r h0 = r h d = r h hd = 0, r h h = 1, r d0 = r d d = r d hd = 0 and r d h = 1. These values are taken fromestimates provided in [19], which are justified therein. Timesare in units of minutes, and rates have units of min − .FIG. 3. Oscillations in Hes/Her protein numbers for the sim-plified 2-cell gene regulatory system with Delta-Notch cou-pling. The system parameters are identical to Fig. 2, but here r h hd = 0 . r h h = 0 .
1, i.e. the coupling between the two cellshas been increased (see Supplement Section S1). As a result,the oscillations in the individual-based system are more clear,periodic and synchronised. This is further demonstrated bythe corresponding power spectra in Fig. 4(c2). Again, thedeterministic trajectory poorly reflects the dynamics of theindividual-based system.
B. Delta-Notch signalling mitigates inhomogeneityand promotes a robust and reliable segmentationclock in noisy oscillators
Having introduced the concept of noise-induced cyclesand the mathematical tools that we will use to analysethem, we now turn our attention to the effect that in-creasing the Delta-Notch signalling strength has on thesenoisy oscillations.It has been observed experimentally [14] that muta-tions in the delta gene give rise to defects in the for-mation of somites. This has been attributed to a de-creased coupling between the cells arising from the mu-tation which, due to slight inhomogeneities between cellsand the stochastic nature of the cellular cycles, leads tothe genetic oscillations in neighbouring cells becomingasynchronous [16, 18]. Furthermore, it has been shownthat encumbered Delta-Notch signalling (i.e. reducedsignalling strength) can give rise to greater disparitiesbetween the oscillations in cells which would be synchro-nised if signalling were not impaired [52]. In this section,we reproduce and study this effect with our model andtheoretical approach (presented in the Supplemental ma-terial in Section S2), thus verifying the necessity of Delta-Notch signalling for the somite segmentation clock. Ajustification for our mathematical definition of ‘couplingstrength’ is given in Supplement Section S1.Consider a two-cell system in which each cell hasslightly different internal parameters (e.g. transcriptionaldelay time) such that the typical cycle time varies be-tween the cells when they are uncoupled. We evaluate,using our theoretical approach, the response of the peakfrequencies, inter-cell phase lag, amplification and coher-ence of genetic oscillations in this inhomogeneous two-cell system for various degrees of Delta-Notch couplingstrength and thus show that the quality of the oscillations(and therefore the segmentation clock) can improve whenDelta-Notch signalling is enhanced.We use two different theoretical approaches in the fol-lowing sections, each of which is valid for different, butcomplementary, parameter sets. (1) In the regime wherethe deterministic (noiseless) equations approach a fixedpoint, we evaluate the power spectrum of fluctuationsof the emergent noise-induced cycles using the so-calledlinear-noise approximation (discussed in more detail Sec-tion S2 of the Supplemental Material). This analysis re-lies upon the deterministic dynamics tending towards afixed point, and approximates the stochastic equationsas linear in the vicinity of this fixed point. The accuracyof the approximation is tested against individual-basedsimulations of the full non-linear stochastic model in Fig.4. (2) When the fixed point of the deterministic systembecomes unstable, we use a deterministic linear stabilityanalysis (LSA) to find the dominant oscillatory frequencyof the cycles and the inter-cell phase lag. The LSA pro-vides no way of finding the amplitude or the coherence ofthe cycles however – it is only possible for us to evaluatethese when method (1) is valid.It has been suggested that a sufficiently strong non-linearity (the so-called ‘cooperativity’) is required for reg-ular sustained cycles [38]. Noting that our linear analysisagrees with simulation results (see Fig. 4), we stress thatthe precise form of the non-linearity is not directly im-portant for the emergence of noise-induced cycles. Theirproperties are well captured by the linearised dynamics.Non-linearities will however affect the location and na-ture of deterministic fixed points, and the coefficients inthe linearised equations near these fixed points.
1. Delta-notch signalling synchronises noisy geneticoscillators and reduces their phase difference
Firstly, we find that Delta-Notch signalling can havethe effect of synchronising the dominant oscillatory fre-quencies of two cells with differing internal parameters.This is demonstrated in Fig. 4, which depicts the powerspectra of the stochastic fluctuations in two such cells forvarious degrees of coupling strength. One observes thatas the coupling strength is increased, the peaks for ei-ther cell, which are separated when there is no coupling,are both drawn towards a common frequency, indicat-ing synchronisation. The degree of synchronisation variessmoothly with the variation of the coupling strength, asshown in Fig. 5(a); in this figure, the dependence ofthe dominant frequencies on coupling strength is shownin more detail. The common frequency which is con-verged upon at large strengths of the Delta-Notch cou- pling agrees with that predicted by linear stability analy-sis (LSA) (see Supplement Section S3 and Section III B 2for further details).Secondly, we find that the peaks of the power spectrain either cell converge to a frequency which reduces thephase lag between the oscillations in the two cells (seeFig. 4). So, not only can Delta-Notch signalling act toalign the frequency of oscillations in neighbouring cells, itcan also encourage the oscillators in either cell to be morealigned in phase. The smooth decrease of the phase lagwith increasing coupling strength is shown in Fig. 5(b).In a similar way to the peak frequency, the phase lagbetween cells agrees with that predicted by LSA whenthe coupling is large.Both of these factors, a shared oscillatory frequencyand a minimal phase lag, are important for the properfunctioning of a cellular clock. The changes in the peakfrequencies and the phase lag that result from an in-crease in coupling strength correspond to quite a notice-able change in the quality of the oscillations themselves.Figs. 2 and 3 are evaluated for the same sets of param-eters as Figs. 4(a) and 4(c) respectively. One parameterset is without coupling between the cells (Fig. 2) andone is with cell-to-cell coupling (Fig. 3). In the formercase, the oscillations in either cell are somewhat aperi-odic and there is no noticeable synchronisation betweenthe cells. However, in Fig. 3, the highs and lows of thecycles in either cell are more inclined to align – this is as-sociated with reduced phase lag and synchronised peakfrequencies. − π πφ ( ω ) (a1) (b1) (c1) .
07 0 . . P ( ω ) (a2) Cell 1Cell 2 .
07 0 . . ω (b2) .
07 0 . . (c2) FIG. 4. Synchronisation of stochastic oscillations of hes/her expression in the 2-cell system as coupling strength is increased.The system parameters are as in Fig. 2 but with (a) r hhd = 0, (b) r hhd = 0 . r hhd = 0 .
9, subject to the constraint r hh + r hhd = 1. That is, the Delta-Notch coupling strength increases from (a) to (c). Panels (a1), (b1) and (c1) depict the phaselag φ ( p ) h , ( ω ) between the oscillations of protein numbers in the two cells as a function of frequency ω . Panels (a2), (b2) and(c2) show the associated Fourier power spectra P ( p ) hj ( ω ) for both cells. In all panels, simulation results are represented bycoloured markers whereas theory results are shown as black lines. The theory lines are produced using the analysis presentedin Sections S2B and S2C of the Supplemental Material. Simulation results are averaged over 100 realisations of the system. Inpanel (a2), there is zero coupling and the peaks of the power spectra are separate, indicating different frequencies of oscillationin the two cells and a lack of synchronisation. In panel (a1), the phase lag between the two cells is random since the two cellsoscillate independently. One observes that as the coupling strength is increased, the cells converge on a common frequency (i.e.they synchronise) and that this common frequency is one which minimises the phase lag between the cells. The power spectrain panels (a2) and (c2) correspond to the time series in Figs. 2(a) and 3(a) respectively.
2. Delta-Notch coupling increases the amplitude andcoherence of noisy oscillations
We asserted previously that an important character-istic of an effective cellular clock is a well-defined time-period – if many frequencies contribute significantly tothe oscillations, then it is more difficult to identify anoverall phase for the clock. We also asserted the neces-sity of the cycles to have a significant amplitude. Bothof these factors contribute to the clarity of the ‘signal’of the oscillations that constitute the cellular clock. Wenote that amplification and coherence are properties ofthe cycles in individual cells whereas synchronisation andphase lag are comparative measures of the oscillations indifferent cells. Despite this, amplification and coherenceare indeed affected by Delta-Notch coupling too.We find that as the Delta-Notch coupling strength isincreased, the amplitude of the oscillations in both cellsfirst decreases slightly then increases, as shown in Fig.5(c). A similar result was also found in previous ex-perimental and theoretical works [17, 35]. However, itis important to note the following caveat on the resultspresented in Fig. 5(c):The blue and red lines shown in Fig. 5 were produced using the theoretical approach based on the linear-noiseapproximation (see Sec. S2 of the Supplemental mate-rial). As such, we observe that the calculated amplitudeof the oscillations diverges when the coupling strength isincreased sufficiently. This singularity is a consequence ofour theoretical approximation; it corresponds to the on-set of an instability for the fixed point of the determin-istic (noiseless) equations and the emergence of a limitcycle [62, 63]. The nonlinearities in the model then be-come more relevant and curtail the amplitude of the os-cillations; this is not captured by the linear theoreticalapproach. The point at which this instability occurs ispredicted by deterministic LSA (see Supplement SectionS3) and is indicated by dashed vertical lines in Fig. 5.It is at this point that the stochastic theory, which isvalid only when the deterministic system approaches afixed point, becomes inaccurate. Instead, the LSA be-comes the more accurate tool for identifying the peakfrequency of oscillation and the inter-cell phase lag. Thetwo analytical methods complement each other in thissense – we can use both approaches to continuously anal-yse the dominant frequency and phase-lag over the onsetof the deterministic instability. Unfortunately, the LSAprovides no means of calculating the amplification or thecoherence of the cycles – we are only able to accuratelypredict these quantities when the deterministic dynam-ics approaches a fixed point (i.e. before the onset of thedeterministic instability).With this caveat in mind, one can nevertheless see fromFig. 4 that the linear stochastic theory still agrees withsimulation results close to this transition. That is, itremains accurate over a wide enough range of couplingstrengths to faithfully capture an increase in the ampli-fication with coupling strength, which occurs before theonset of the instability.We find also that as the coupling strength is increased,the power spectrum of fluctuations becomes sharplypeaked (as can be seen in Fig. 4(c)) at a characteris-tic frequency, i.e., the cycles become more coherent (seeFig. 5(d)). The location of this peak in the power spec-trum corresponds to the frequency to which the cyclesin the two cells converge (as discussed in the previoussection).We conceptualise this increase in amplification and co-herence as a consequence of a kind of ‘resonant amplifica-tion’. 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 theinter-cell coupling strength is increased, the frequenciesare aligned and the phase lag between them is reduced,the cycles begin to constructively interfere at a charac-teristic frequency.Interestingly, as the coupling strength is increased from zero, the amplitude of the oscillations in either cell ini-tially decreases (Fig. 5(c)), as does the peakedness ofthe power spectrum (Fig. 5(d)). We attribute this to thefact that, for low coupling, the oscillations in either cellare not adequately synchronised for their interference tobe constructive. So as the coupling strength is increasedinitially, the ‘interference’ between the two cells has a de-structive effect. It is only when the phase lag is reducedand the frequencies are aligned sufficiently (as a resultof a further increase of the coupling) that the collectiveamplitude of oscillations increases.The effects of increased amplification and coherenceare evident in the time-series of the noisy cellular os-cillations in the two-cell system shown in Figs. 2 and3 (which are evaluated with and without inter-cellularcoupling respectively). In Fig. 2, there is no easily dis-cernible periodic nature to the cycles in either cell. Thisis in contrast to the cycles in Fig. 3 where the highs andlows have more consistent temporal separations. Thisgreater clarity in the oscillatory frequency is associatedwith the the increase in the sharpness of the peaks ofthe power spectra between Figs. 4(a2) and 4(c2) (i.e.an increase in coherence). Additionally, it can be seenthat there are fewer pronounced highs and lows, on thewhole, in the uncoupled system than in the coupled sys-tem; in Fig. 2 highs and lows are sporadically interruptedby stints of somewhat suppressed fluctuations about thefixed point. This in turn contributes to the lower overallamplification of the uncoupled system in comparison tothat of the coupled system.
FIG. 5. Synchronisation (a), phase lag (b), amplification (c) and coherence (d) of oscillations of Hes/Her protein in theinhomogeneous two-cell system versus coupling strength produced using the linear theory detailed in Section S2 of the Sup-plemental Material. We use ∆ ω = 0 .
01 in our definition of coherence (see Eq. (5)). The system parameters are the same asin Fig. 2 but r hhd and r hh are varied subject to the constraint r hh + r hhd = 1 (see Supplement Section S1 for the definition ofthese parameters). Panel (a) demonstrates that the peak frequency of oscillation ω ( j )max in either cell (where j = 1 , r hhd is increased. This limiting value agrees with thatpredicted by linear stability analysis in the deterministic system (see Supplement Section S3). Panel (b) shows that as the twocells synchronise, the (rescaled) phase lag φ , ( ω ( j )max ) /π between the two cells also decreases. Again, the value of the phaselag agrees with linear stability analysis for strong coupling. Panel (c) demonstrates how the oscillations in the two cells areinitially dampened and subsequently amplified as the coupling strength is increased from zero. The theory becomes invalid asthe deterministic fixed point of the system becomes unstable but is accurate close up to this point (see main text). Panel (d)shows how the coherence of the power spectra initially decreases and then increases as the coupling strength is increased, in asimilar way to the amplification shown in panel (c). Values of the red and blue lines close to and to the right of the verticaldotted line at r hhd = 0 .
985 in all panels are not necessarily accurate; in this regime the system is close to or beyond the onsetof the deterministic instability indicated by the vertical dashed line (emergence of a limit cycle). This results in corrections tothe Fourier spectra which are not accounted for in our linear theory (see Supplement). It is at this point that the LSA, whichis a deterministic analysis, becomes useful for identifying the peak frequency and phase lag. For low coupling strength (beforethe instability), the LSA is inaccurate because it does not take into account the effects of the noise, highlighting the need forthe stochastic theory.
C. Transcriptional/translational delays can lead toout-of-phase oscillations despite Delta-Notchcoupling
We demonstrated in the previous sections that twocells with slightly disparate oscillatory frequencies couldsynchronise when coupled via Delta-Notch signalling. Asthe strength of the Delta-Notch coupling is increased, acommon frequency is converged upon and the phase lag between the two cells is reduced. Although this is char-acteristic of many model parameter sets, it is not alwaysthe case. As was noted also in [19] and [31] (in the purelydeterministic setting), neighbouring cells can be coupledin such a way that they oscillate in anti-phase with oneanother. This type of behaviour is facilitated by the de-lays associated with translation and/or transcription.Fig. 6 demonstrates that for certain parameter sets,as the coupling is increased, oscillations in two cells can0converge on a common frequency but the phase lag be-tween the two cells can approach values closer to φ = π than φ = 0. That is, the cells tend towards oscillating in anti-phase with one another. Clearly, this is suboptimalif these cellular oscillations are to be used as a segmen-tation clock. FIG. 6. Peak frequencies ω ( j )max and rescaled phase lag φ , ( ω ( j )max ) /π in the two-cell system (where j = 1 , r hhd and r hh are varied subject to the constraint r hh + r hhd = 1 and τ ( m ) d = 20. Panel (a) demonstratesthat as the coupling strength r hhd is increased, the peak frequency of oscillations approaches a common value, as in Fig. 5.However, panel (b) shows that as the two cells synchronise, the phase lag between the two cells increases, tending towards π instead of 0. To understand why this should happen, we monitor the dominant frequency and the associated phase differencebetween the two cells as the transcriptional/translational delays are varied (see Fig. 7). For the purposes of thisanalysis, the two cells are taken to be identical. In this case, the two cells are guaranteed to share a peak oscillatoryfrequency ω (1)max = ω (2)max = ω max . We find that whether the cells oscillate in or out-of phase is determined by aninterplay between the delays and the dominant frequency of oscillation. FIG. 7. Half the time period associated with dominant frequency of oscillation in the 2-cell system
T / π/ω max (a) andthe corresponding phase lag between the cells (b). The system parameters are the same as in Fig. 3 but τ ( m ) h = 20 in bothcells and τ ( p ) d and τ ( m ) d are varied. The typical half-period varies around an average value of ∼ T / τ ( p ) d and τ ( m ) d , it typically remains within ∼
10% of this mean value. The phase lag switchesbetween 0 and π along the lines lines τ ( p ) d + τ ( m ) d = τ tot i (where i = 1 , , , . . . ). The values τ tot i at which the switch in phaseoccurs are separated by regular intervals such that τ tot i − τ tot i − ≈ One observes that the phase lag between the two cells switches between 0 and π (and vice versa) when the to-1tal delta delay time τ tot = τ ( m ) d + τ ( p ) d reaches certainvalues τ tot i , where i = 1 , , , . . . (we notice a similar ef-fect when other pairs of delays are varied). In Figure 7,the times at which the switches occur are separated bya regular time interval such that τ tot i − τ tot i − ≈ T / π/ω max , which varies withinaround ∼
10% of the mean value ∼ π/ω max to the total delay time, the ‘signal’ from one cell to itsneighbour is delayed by half a cycle. If the cells would os-cillate in phase without this additional delay, it stands toreason that they would oscillate in anti-phase given theadditional delay – there would be no difference from thepoint of view of either cell (assuming that the frequencyof oscillation does not vary greatly with the changing de-lay times). A similar argument holds true for the additionof 2 π/ω max to τ tot – in this case the phase difference be-tween the two cells ought to be unchanged. One caveatto this reasoning is that the frequency of the oscillationscannot be changed too drastically by the variation in timedelay.We conclude that while the time delays associated withtranslation and transcription are crucial for the persis-tence of the cycles which constitute the cell clock, theyare somewhat of a double-edged sword. Depending onthe interaction between the delays and the internal clocksof each cell, the cells may oscillate in anti-phase with oneanother. As a result, the precise nature of the transcrip-tional/translational delays is of great importance withregards to the proper synchronisation of the segmenta-tion clock. D. Transcriptional/translational delays can disruptglobal oscillations in chains of cells and give rise towaves of gene expression
In this section, we discuss how the preceding analysiscan be extended to a chain of Delta-Notch-coupled cells.We demonstrate that synchronised noisy oscillations inthe two-cell system can correspond to global oscillationsin a chain. We also explore behaviours other than globaloscillations which can occur as a result of delays; namely,the emergence of noise-induced waves.A distinction between oscillations of the determinis-tic trajectory and purely noise-induced oscillations wasmade in Section III A. In a similar way, one finds thatwaves can manifest in the individual-based system whenthey do not in the deterministic system. Such ‘stochas-tic waves’ or ‘quasi-waves’ have been found previouslyin theoretical models of individual-based systems withlong-range interaction [64]. Here however, the stochas-tic waves arise due to a combination of the non-localdependence of the reaction rates (due to Delta-Notchsignalling) and the transcriptional/translational delays. Conversely, waves of gene expression have been studiedpreviously in chains of coupled genetic oscillators [31] butthis was done in the context of deterministic equationswhich ignored intrinsic noise.We mention the emergence of waves here not so muchas an explanation for the travelling waves which are seenin the PSM (these are most likely due to a variationin translatonal/transcriptional delay along the anterior-posterior axis [35, 52]), but as an illustration of the differ-ent kinds of undesirable behaviour which can arise whencells oscillate out of phase with one another. As such,the results presented in section are exploratory and notnecessarily an attempt to recreate any (as yet) observedphenomena.We find that for sets of parameters where one wouldobserve oscillations in anti-phase in the two-cell system,one finds waves of gene expression in an extended chain ofcells. For parameter sets where the cycles in the two-cellsystem oscillate in unison, one observes global in-phaseoscillations in gene expression.For an extended chain of cells, we define the powerspectrum P αk ( ω ) = (cid:104)| ˆ˜ δ αk ( ω ) | (cid:105) , (6)where we have used the discrete Fourier transformwith respect to the cell number j defined by ˜ f k = √ L (cid:80) j e ijk f j , where L is the number of cells in the chain.Details of the calculation of this power spectrum aregiven in the Supplement (Section S2 B).In a chain of coupled cells, global in-phase oscillationsare characterised by a peak in the power spectrum P αk ( ω )at spatial wavenumber k = 0 and non-zero oscillatoryfrequency ω . Such a power spectrum is shown in Fig.8(a), and an example of the corresponding behaviour ina chain of cells is demonstrated in Fig. 8(b), where oneobserves that the peaks and troughs in one cell tend toalign with those in the neighbouring cells. That the cellsare indeed in phase with one another is verified by thephase lag (see inset of Fig. 8(a)) which is equal to zero,regardless of cell separation.Stochastic waves, on the other hand, are characterisedby a peak in the Fourier power spectrum at non-zero val-ues of both the spatial wavenumber k and the temporalfrequency ω . An example of such a power spectrum isgiven in Fig. 9(a). To validate the claim that this peakin the power spectrum is indicative of travelling waves,we note that the phase difference between cells varies lin-early with cell separation, as is shown in the inset of Fig.9(a). We stress that the coupling between the cells hereis biased in one direction, which breaks the symmetry ofthe system, allowing waves to travel. Travelling wavesof gene expression have been observed in experimentalsystems other than the PSM [65, 66].For symmetric coupling however, one instead observesstanding waves of gene expression, where alternate cellsoscillate in antiphase, as shown in Fig. 10. In this partic-ular case, the phase lag between any pair of adjacent cells2(at the peak frequency) is π , as can be seen in the inset ofFig. 10(a). This is rather reminiscent of the on-off che-querboard patterns associated with neural differentiationand lateral inhibition [67].Examples of travelling and standing stochastic wavesin a chain of coupled cells are shown in Figs. 9(b) and10(b) respectively. There is a clear qualitative distinction between the two. In Fig. 9(b), peaks and troughs inone cell gradually travel in the positive j direction astime goes by, indicating a travelling wave. In Fig. 10(b)however, the peaks in one cell tend to line up with thetroughs of the neighbouring cells and visa versa and thereis no clear direction of travel. .
09 0 .
10 0 .
11 0 .
12 0 . Temporal frequency ω π π π S pa t i a l w a v enu m be r k (a) P o w e r s pe c t r u m P ( p ) h k ( ω ) | j − j |− π πφ j,j Time t P o s i t i on j (b) P r o t e i nnu m be r n ( p ) h j ( t ) FIG. 8. Global oscillations in Hes/Her protein numbers. Panel (a): The power spectrum of fluctuations P ( p ) hk ( ω ) and thecorresponding phase lag between cells at peak frequency φ ( p ) hj,j (cid:48) ( ω max ) as a function of separation (inset). Global oscillationsare indicated by the fact that the peak of the power spectrum is located at k = 0 and non-zero ω and also by the phase lagbetween cells being zero. Panel (b): An example of noise-induced global oscillations. Data is from a simulation of the stochasticindividual-based model. The peaks and troughs in hes/her expression have a tendency to align, giving rise to vertical stripedstructures in the figure. The coupling between cells is symmetric ( d (+) = d ( − ) = 1). The remaining system parameters are thesame as in Fig. 3 with the exception that τ ( m ) h = 20 in all cells. .
09 0 .
10 0 .
11 0 .
12 0 . Temporal frequency ω π π π S pa t i a l w a v enu m be r k (a) P o w e r s pe c t r u m P ( p ) h k ( ω ) j − j π π πφ j,j Time t P o s i t i on j (b) P r o t e i nnu m be r n ( p ) h j ( t ) FIG. 9. Travelling stochastic waves. Panel (a): The power spectrum of fluctuations P ( p ) hk ( ω ) and the corresponding phaselag between cells at peak frequency φ ( p ) hj,j (cid:48) ( ω max ) as a function of cell separation (inset). The presence of travelling waves isindicated by a sharp peak in the power spectrum at a non-zero value of ω and a value of k that is neither equal to 0 nor π .The presence of travelling waves is further evidenced by a phase lag which is linearly increasing with cell separation. We notethat we have relaxed the condition that φ ( p ) hj,j (cid:48) be in the range [ − π, π ) in order for this trend to be apparent. Panel (b): Anexample of noise-induced travelling waves; this is from simulations of the stochastic individual-based model. When a peak ora trough occurs in one cell, it has a tendency to move upwards to the neighbouring cell as time progresses, giving rise to thediagonal structures in the figure. The coupling between cells is asymmetric ( d (+) = 2, d ( − ) = 0), which breaks the directionalsymmetry of the system, allowing waves to propagate. The remaining system parameters are the same as in Fig. 2 with theexception that r hh = 0 . r hhd = 0 . τ ( m ) h = 20 in all cells and τ ( m ) d = 35. .
09 0 .
10 0 .
11 0 .
12 0 . Temporal frequency ω π π π S pa t i a l w a v enu m be r k (a) P o w e r s pe c t r u m P ( p ) h k ( ω ) | j − j | πφ j,j Time t P o s i t i on j (b) P r o t e i nnu m be r n ( p ) h j ( t ) FIG. 10. Neighbouring cells oscillating in anti-phase (standing waves). Panel (a): The power spectrum of fluctuations P ( p ) hk ( ω )and the corresponding phase lag between cells at peak frequency φ ( p ) hj,j (cid:48) ( ω max ) as a function of separation (inset). Neighbouringcells tend to oscillate in anti-phase with one another, as demonstrated by the phase lag (inset) and also by the location ofthe peak of the power spectrum at k = π and non-zero ω . Panel (b): An example of noise-induced standing waves seen insimulations of the stochastic individual-based model. The peaks in one cell tend to align with the troughs in the neighbouringcell, giving rise to chequerboard-type patterns in the figure. The coupling between cells is symmetric ( d (+) = d ( − ) = 1). Theremaining system parameters are the same as in Fig. 2 with the exception that r hh = 0 . r hhd = 0 . τ ( m ) h = 20 in all cellsand τ ( m ) d = 20. IV. SUMMARY AND DISCUSSION
The process of somite segmentation poses a complexmany-faceted problem for theorists and experimentalists alike and remains an area of active inquiry. Much is leftto be discovered about the precise nature of the role ofeach of the genes involved with the somite segmentation4clock and their interactions with external signalling fac-tors in the embryo.In this paper, we aim to have provided some insightinto the role that Delta-Notch signalling plays in not onlyaligning the frequencies of oscillation of cyclic gene ex-pression in neighbouring cells but also in reducing phaselag, in improving the coherence of oscillations and in in-creasing their amplitude, so as to produce a robust andreliable segmentation clock. We also explored the rolethat intrinsic noise plays in the system; counter to in-tuition, it can actually promote persistent cycles, ratherthan obscure them. Further, we discussed how the de-lays involved in the transcription and translation pro-cesses can act to promote oscillations but can also resultin neighbouring cells oscillating out-of-phase with one an-other. We examined how this resulted from an interplaybetween the dominant frequency of oscillation in the cellswith the aggregate time-delay. We went on to show howasynchronous behaviour in a two-cell system correspondsto waves of gene expression in a chain of cells.In a recent work [52], the gene expression noise in thePSM of the zebrafish was analysed. This was done by us-ing smFISH microscopy techniques to count the numbersof discrete RNA molecules in individual cells. The statis-tical discrepancies between the gene expression in sets ofcells which were supposedly synchronised was then eval-uated (using well-known techniques [27, 28]) and termed‘expression noise’. It was found that this expression noiseincreased 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 in-creased Delta-Notch signalling strength increases the de-gree of synchronisation and the coherence of noisy oscil-lations, which in turn reduces the discrepancy betweenthe cycles in coupled cells. This would appear to be verymuch in keeping with the aforementioned experimental findings.On a more general note, intrinsic noise is often assumedto be a destabilising influence on cycles and to be thesource of the discrepancy between the expression in twootherwise-equal cells. But, because of the complexity ofthe gene regulatory network and the nature of the cou-pling between cells in the PSM, the intrinsic noise canactually give rise to correlated sustained oscillations inneighbouring cells– a behaviour one might normally as-sociate with an extrinsic influence. This rather blurs theline between the what might be considered the signaturesof intrinsic, extrinsic and ‘expression’ noise in experimen-tal data. As a result the utmost care must be taken toidentify sources of correlation between cells, other thancommon external influence, if one is to truly discern thefingerprints of intrinsic noise in the data from those ofextrinsic stochasticity.
AUTHOR CONTRIBUTIONS
JWB designed the study, contributed to discussionsguiding the work, carried out mathematical calculations,performed simulations and wrote the manuscript. TGdesigned the study, contributed to discussions guidingthe work and wrote the manuscript.
FUNDING STATEMENT
JWB thanks the Engineering and Physical SciencesResearch Council (EPSRC) for funding (PhD stu-dentship, EP/N509565/1). TG acknowledges partial fi-nancial support from the Maria de Maeztu Program forUnits of Excellence in R&D (MDM-2017-0711). We de-clare no competing interests. [1] Olivier Pourqui´e, “The segmentation clock: convertingembryonic time into spatial pattern,” Science , 328–330 (2003).[2] Andrew C. Oates, Luis G. Morelli, and Sa´ul Ares, “Pat-terning embryos with oscillations: structure, functionand dynamics of the vertebrate segmentation clock,” De-velopment , 625–639 (2012).[3] Yumiko Saga and Hiroyuki Takeda, “The making of thesomite: molecular events in vertebrate segmentation,”Nat. Rev. Genet. , 835 (2001).[4] John Cooke and Erik Christopher Zeeman, “A clock andwavefront model for control of the number of repeatedstructures during animal morphogenesis,” J. Theor. Biol. , 455–476 (1976).[5] Yasumasa Bessho, Ryoichi Sakata, Suguru Komatsu, Ko-hei Shiota, Shuichi Yamada, and Ryoichiro Kageyama,“Dynamic expression and essential functions of hes7 insomite segmentation,” Genes Dev. , 2642–2647 (2001).[6] Clarissa A Henry, Michael K Urban, Kariena K Dill,John P Merlie, Michelle F Page, Charles B Kimmel, and Sharon L Amacher, “Two linked hairy/enhancer of split-related zebrafish genes, her1 and her7, function togetherto refine alternating somite boundaries,” Development , 3693–3704 (2002).[7] Andrew C Oates and Robert K Ho, “Hairy/e (spl)-related (her) genes are central components of the seg-mentation oscillator and display redundancy with thedelta/notch signaling pathway in the formation of ante-rior segmental boundaries in the zebrafish,” Development , 2929–2946 (2002).[8] Dirk Sieger, Bastian Ackermann, Christoph Winkler, Di-ethard Tautz, and Martin Gajewski, “her1 and her13. 2are jointly required for somitic border specification alongthe entire axis of the fish embryo,” Dev. Biol. , 242–251 (2006).[9] Ronald A Conlon, Andrew G Reaume, and JanetRossant, “Notch1 is required for the coordinate segmen-tation of somites,” Development , 1533–1545 (1995).[10] Martin Hrabˇe de Angelis, Joseph Mclntyre II, and AchimGossler, “Maintenance of somite borders in mice requires the delta homologue dll1,” Nature , 717 (1997).[11] Kenro Kusumi, Eileen Sun, Anne W Kerrebrock, Rod-erick T Bronson, Dow-Chung Chi, Monique Bulotsky,Jessica B Spencer, Bruce W Birren, Wayne N Frankel,and Eric S Lander, “The mouse pudgy mutation dis-rupts delta homologue dll3 and initiation of early somiteboundaries,” Nat. Genet. , 274 (1998).[12] Alexander Aulehla, Christian Wehrle, Beate Brand-Saberi, Rolf Kemler, Achim Gossler, Benoit Kanzler,and Bernhard G Herrmann, “Wnt3a plays a major role inthe segmentation clock controlling somitogenesis,” Dev.Cell , 395 – 406 (2003).[13] Julien Dubrulle, Michael J. McGrew, and OlivierPourqui, “Fgf signaling controls somite boundary po-sition and regulates segmentation clock control of spa-tiotemporal hox gene activation,” Cell , 219 – 232(2001).[14] C. Jouve, I. Palmeirim, D. Henrique, J. Beckers,A. Gossler, D. Ish-Horowicz, and O. Pourquie, “Notchsignalling is required for cyclic expression of the hairy-like gene hes1 in the presomitic mesoderm,” Development , 1421–1429 (2000).[15] Scott A. Holley, Robert Geisler, and Christiane Nsslein-Volhard, “Control of her1 expression during zebrafishsomitogenesis by a delta-dependent oscillator and an in-dependent wave-front activity,” Genes Dev. , 1678–1690 (2000).[16] Yun-Jin Jiang, Birgit L Aerne, Lucy Smithers, CatherineHaddon, David Ish-Horowicz, and Julian Lewis, “Notchsignalling and the synchronization of the somite segmen-tation clock,” Nature , 475 (2000).[17] Ertu˘grul M ¨Ozbudak and Julian Lewis, “Notch signallingsynchronizes the zebrafish segmentation clock but is notneeded to create somite boundaries,” PLOS Genetics ,1–11 (2008).[18] Emilie A. Delaune, Paul Franois, Nathan P. Shih, andSharon L. Amacher, “Single-cell-resolution imaging ofthe impact of notch signaling and mitosis on segmenta-tion clock dynamics,” Developmental Cell , 995 – 1005(2012).[19] Julian Lewis, “Autoinhibition with transcriptional delay:a simple mechanism for the zebrafish somitogenesis oscil-lator,” Curr. Biol. , 1398–1408 (2003).[20] Ertu˘grul M ¨Ozbudak, Mukund Thattai, Iren Kurtser,Alan D Grossman, and Alexander Van Oudenaarden,“Regulation of noise in the expression of a single gene,”Nature genetics , 69 (2002).[21] DW Austin, MS Allen, JM McCollum, RD Dar, JR Wil-gus, GS Sayler, NF Samatova, CD Cox, and ML Simp-son, “Gene network shaping of inherent noise spectra,”Nature , 608 (2006).[22] Hiromi Shimojo, Akihiro Isomura, Toshiyuki Ohtsuka,Hiroshi Kori, Hitoshi Miyachi, and Ryoichiro Kageyama,“Oscillatory control of delta-like1 in cell interactions reg-ulates dynamic gene expression and tissue morphogene-sis,” Genes Dev. , 102–116 (2016).[23] Nick E Phillips, Cerys S Manning, Tom Pettini, Veron-ica Biga, Elli Marinopoulou, Peter Stanley, James Boyd,James Bagnall, Pawel Paszek, David G Spiller, et al. ,“Stochasticity in the mir-9/hes1 oscillatory network canaccount for clonal heterogeneity in the timing of differ-entiation,” eLife , e16118 (2016). [24] Nick E Phillips, Cerys Manning, Nancy Papalopulu, andMagnus Rattray, “Identifying stochastic oscillations insingle-cell live imaging time series using gaussian pro-cesses,” PLOS Comput. Biol. , e1005479 (2017).[25] Hiromi Hirata, Shigeki Yoshiura, Toshiyuki Ohtsuka, Ya-sumasa Bessho, Takahiro Harada, Kenichi Yoshikawa,and Ryoichiro Kageyama, “Oscillatory expression of thebhlh factor hes1 regulated by a negative feedback loop,”Science , 840–843 (2002).[26] Hiromi Shimojo, Toshiyuki Ohtsuka, and RyoichiroKageyama, “Oscillations in notch signaling regulatemaintenance of neural progenitors,” Neuron , 52 – 64(2008).[27] Michael B. Elowitz, Arnold J. Levine, Eric D. Siggia, andPeter S. Swain, “Stochastic gene expression in a singlecell,” Science , 1183–1186 (2002).[28] Peter S. Swain, Michael B. Elowitz, and Eric D. Siggia,“Intrinsic and extrinsic contributions to stochasticity ingene expression,” Proc. Natl. Acad. Sci. , 12795–12800(2002).[29] Peter B. Jensen, Lykke Pedersen, Sandeep Krishna, andMogens H. Jensen, “A wnt oscillator model for somito-genesis,” Biophys. J. , 943 – 950 (2010).[30] Albert Goldbeter and Olivier Pourqui, “Modeling thesegmentation clock as a network of coupled oscillationsin the notch, wnt and fgf signaling pathways,” J. Theor.Biol. , 574 – 585 (2008), in Memory of Reinhart Hein-rich.[31] Hiroshi Momiji and Nicholas A. M. Monk, “Oscillatorynotch-pathway activity in a delay model of neuronal dif-ferentiation,” Phys. Rev. E , 021930 (2009).[32] Olivier Cinquin, “Is the somitogenesis clock really cell-autonomous? a coupled-oscillator model of segmenta-tion,” J. Theor. Biol. , 459 – 468 (2003).[33] Luis G Morelli, Sa´ul Ares, Leah Herrgen, ChristianSchr¨oter, Frank J¨ulicher, and Andrew C Oates, “Delayedcoupling theory of vertebrate segmentation,” HFSP J. ,55–66 (2009).[34] Kazuki Horikawa, Kana Ishimatsu, Eiichi Yoshimoto,Shigeru Kondo, and Hiroyuki Takeda, “Noise-resistantand synchronized oscillation of the segmentation clock,”Nature , 719 (2006).[35] Ahmet Ay, Jack Holland, Adriana Sperlea, Gnanapack-iam Sheela Devakanmalai, Stephan Knierer, SebastianSangervasi, Angel Stevenson, and Ertu˘grul M. ¨Ozbudak,“Spatial gradients of protein-level time delays set thepace of the traveling segmentation clock waves,” Devel-opment , 4158–4167 (2014).[36] Ruth E. Baker, Santiago Schnell, and Philip K. Maini,“Mathematical models for somite formation,” in Multi-scale Modeling of Developmental Systems , Current Top-ics in Developmental Biology, Vol. 81 (Academic Press,2008) pp. 183 – 203.[37] Tobias Galla, “Intrinsic fluctuations in stochastic delaysystems: Theoretical description and application to asimple model of gene regulation,” Physical Review E ,021909 (2009).[38] Manuel Barrio, Kevin Burrage, Andr´e Leier, andTianhai Tian, “Oscillatory regulation of hes1: discretestochastic delay modelling and simulation,” PLoS com-putational biology , e117 (2006).[39] Robert Schlicht and Gerhard Winkler, “A delay stochas-tic process with applications in molecular biology,” J.Math. Biol. , 613–648 (2008). [40] Thomas E Turner, Santiago Schnell, and Kevin Burrage,“Stochastic approaches for modelling in vivo reactions,”Comput. Biol. Chem. , 165–178 (2004).[41] Mukund Thattai and Alexander Van Oudenaarden, “In-trinsic noise in gene regulatory networks,” Proc. Natl.Acad. Sci. , 8614–8619 (2001).[42] Joe H. Levine, Yihan Lin, and Michael B. Elowitz,“Functional roles of pulsing in genetic circuits,” Science , 1193–1200 (2013).[43] Albert Goldbeter, Biochemical oscillations and cellularrhythms: the molecular bases of periodic and chaotic be-haviour (Cambridge university press, 1997).[44] Sarah J Bray, “Notch signalling: a simple pathway be-comes complex,” Nat. Rev. Mol. Cell Biol. , 678 (2006).[45] Emma R. Andersson, Rickard Sandberg, and UrbanLendahl, “Notch signaling: simplicity in design, versa-tility in function,” Development , 3593–3612 (2011).[46] Olivier Pourqui, “Vertebrate segmentation: From cyclicgene networks to scoliosis,” Cell , 650 – 663 (2011).[47] C. W. Gardiner, Handbook of stochastic methods forphysics, chemistry and the natural sciences , 4th ed.(Springer, New York, 2009).[48] Hal L Smith,
An introduction to delay differential equa-tions with applications to the life sciences (Springer, NewYork, 2011).[49] Joseph D. Challenger and Alan J. McKane, “Synchro-nization of stochastic oscillators in biochemical systems,”Phys. Rev. E , 012107 (2013).[50] G. Rozhnova, A. Nunes, and A. J. McKane, “Phase lagin epidemics on a network of cities,” Phys. Rev. E ,051912 (2012).[51] David Alonso, Alan J McKane, and Mercedes Pascual,“Stochastic amplification in epidemics,” J. Royal Soc. In-terface , 575–582 (2006).[52] Sevdenur Keskin, Gnanapackiam S. Devakanmalai,Soo Bin Kwon, Ha T. Vu, Qiyuan Hong, Yin Yeng Lee,Mohammad Soltani, Abhyudai Singh, Ahmet Ay, andErtu˘grul M. ¨Ozbudak, “Noise in the vertebrate segmen-tation clock is boosted by time delays but tamed by notchsignaling,” Cell Rep. , 2175 – 2185.e4 (2018).[53] Ryoichiro Kageyama, Yasutaka Niwa, Akihiro Isomura,Aitor Gonz´alez, and Yukiko Harima, “Oscillatory geneexpression and somitogenesis,” Wiley Interdiscip. Rev.Dev. Biol. , 629–641 (2012).[54] Carlos A Lugo and Alan J McKane, “Quasicycles in aspatial predator-prey model,” Phys. Rev. E , 051911(2008).[55] Alan J McKane and Timothy J Newman, “Predator-prey cycles from resonant amplification of demographicstochasticity,” Phys. Rev. Lett. , 218102 (2005).[56] Andrew J Black and Alan J McKane, “Stochastic ampli-fication in an epidemic model with seasonal forcing,” J.Theor. Biol. , 85–94 (2010).[57] Tobias Brett and Tobias Galla, “Stochastic processeswith distributed delays: chemical langevin equationand linear-noise approximation,” Phys. Rev. Lett. ,250601 (2013).[58] Tobias Brett and Tobias Galla, “Gaussian approxima-tions for stochastic systems with delay: chemical langevinequation and application to a brusselator system,” J.Chem. Phys. , 124112 (2014).[59] Joseph W. Baron and Tobias Galla, “Stochastic fluctu-ations and quasipattern formation in reaction-diffusionsystems with anomalous transport,” Phys. Rev. E , 052124 (2019).[60] David F. Anderson, “A modified next reaction methodfor simulating chemical systems with time dependentpropensities and delays,” The Journal of ChemicalPhysics , 214107 (2007).[61] Daniel T Gillespie, “Exact stochastic simulation of cou-pled chemical reactions,” J. Phys. Chem. , 2340–2361(1977).[62] Richard P Boland, Tobias Galla, and Alan J McKane,“Limit cycles, complex floquet multipliers, and intrinsicnoise,” Phys. Rev. E , 051131 (2009).[63] Richard P Boland, Tobias Galla, and Alan J McKane,“How limit cycles and quasi-cycles are related in systemswith intrinsic noise,” J. Stat. Mech. Theory Exp. ,P09001 (2008).[64] Tommaso Biancalani, Tobias Galla, and Alan J McKane,“Stochastic waves in a brusselator model with nonlocalinteraction,” Phys. Rev. E , 026201 (2011).[65] B´en´edicte Wenden, David L. K. Toner, Sarah K. Hodge,Ramon Grima, and Andrew J. Millar, “Spontaneousspatiotemporal waves of gene expression from biologicalclocks in the leaf,” Proc. Natl. Acad. Sci. , 6757–6762(2012).[66] Kazuya Ukai, Koji Inai, Norihito Nakamichi, HirokiAshida, Akiho Yokota, Yusuf Hendrawan, HaruhikoMurase, and Hirokazu Fukuda, “Traveling waves of cir-cadian gene expression in lettuce,” Environ. Control Biol. , 237–246 (2012).[67] Joanne R Collier, Nicholas AM Monk, Philip K Maini,and Julian H Lewis, “Pattern formation by lateral in-hibition with feedback: a mathematical model of delta-notch intercellular signalling,” J. Theor. Biol. , 429–446 (1996).[68] P. ´Erdi and J. T´oth, Mathematical models of chemicalreactions: theory and applications of deterministic andstochastic models (Manchester University Press, Manch-ester, 1989).[69] A. Altland and B. D. Simons,
Condensed matter fieldtheory (Cambridge University Press, Cambridge, 2010).[70] N. G. Van Kampen,
Stochastic processes in physics andchemistry , Vol. 1 (Elsevier, London, 1992).[71] D. T. Gillespie, J. Chem. Phys. , 1716 (2001). ntrinsic noise, Delta-Notch signalling and delayed reactions promotesustained, coherent, synchronised oscillations in the presomiticmesodermSupplemental Material
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 transcrip-tion/translation processes and finite particle numbers) using an individual-based model. In this model, we supposethat 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 delaysbetween 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, M d a d = ⇒ τ ( p ) d M d + P d ,M h a h = ⇒ τ ( p ) h M h + P h , ∅ f d = ⇒ τ ( m ) d M d , ∅ f h = ⇒ τ ( m ) h M h ,P d b d −→ ∅ ,P h b h −→ ∅ ,M d c d −→ ∅ ,M h c h −→ ∅ , (S1)where M d denotes a molecule of delta mRNA and P h denotes a molecule of Hes/Her protein, etc. The single arrows R −→ indicate that the reaction occurs without delay with a per capita rate R . The double arrows R = ⇒ τ indicate adelayed reaction with per capita rate R and delay τ . These equations are to be interpreted in the usual way usingmass action kinetics [68]. For example, the first reaction is triggered with rate a d n ( m ) d , where n ( m ) d is the numberof M d -particles in the system. The effect of such a reaction is realised τ ( p ) d units of time later, and results in theaddition of a P d -particle to the system.The individual-based dynamics summarised by Eq. (S1) can be approximated by a set of stochastic differentialequations (SDEs), which are given in Section S2 C.The model parameters a h , a d , b h , b d and c h , c 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 f αj ( n t ) = k α N (cid:34) r α + r α d 12 (cid:0) d ( − ) φ d j − ,t + d (+) φ d j +1 ,t (cid:1) (cid:0) d ( − ) φ d j − ,t + d (+) φ d j +1 ,t (cid:1) + r α h
11 + φ h j,t + r α hd 12 (cid:0) d ( − ) φ d j − ,t + d (+) φ d j +1 ,t (cid:1) (cid:0) d ( − ) φ d j − ,t + d (+) φ d j +1 ,t (cid:1)
11 + φ h j,t (cid:35) , (S2)where r α , r α d , r α h , r α hd and k α are rate constants, N is the system size, φ αj,t = n ( p ) αj,t n ( p ) α , and where d ( ± ) are positiveconstants such that (cid:0) d (+) + d ( − ) (cid:1) = 1; n ( p ) α are reference values for the protein levels. Superscript or subscriptindices α are placeholders, representing the cases α ∈ { h , d } . This reaction scheme is similar to the one used in [19].The different terms in the function f αj ( n t ) correspond to autonomous activation ( r α ), activation by delta only( r α d ), inhibition by hes/her only ( r α h ) and mixed response to hes/her and delta r α hd . We constrain the parameters r α , r α d , r α h , r α hd always to sum to unity i.e. r α + r α d + r α h + r α hd = 1. The expression inside the square brackets in Eq. (S2)can therefore range between 0 and 1 dynamically, depending on the concentrations of mRNA or protein and on theparameters r α , r α d , r α h , r α hd . As a result, the typical mRNA birth rate is characterised by k α N .In the main text, we examine the change in the behaviour of coupled genetic oscillators as we vary the ‘couplingstrength’. In order to isolate the effect of the Delta-Notch signalling from the hes/her auto-repression, we always set r h0 = r hd = 0 and vary r hhd and r hh such that r hh + r hhd = 1. As r hhd increases, so does the coupling strength but the roleof hes/her remains the same. For zero coupling r hhd = 0. For maximal coupling r hhd = 1. We keep the typical mRNAproduction rates k α constant.2In Sections III A, III B 1, III B 2 and III C, since we only consider a 2-cell model, (cid:0) d ( − ) φ d j − ,t + d (+) φ d j +1 ,t (cid:1) = φ d j (cid:48) ,t where j (cid:48) = 2 if j = 1 and vice versa. S2. QUANTIFICATION OF STOCHASTIC FLUCTUATIONS IN SYSTEMS WITH DELAYS ANDNON-LOCAL REACTION RATES
In this section, we derive the analytical results for the quantification of the stochastic fluctuations about the fixedpoint 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) [47] which approximate the individual-based dynamics in the limit of large system-size N . We then use the linear-noise approximation to obtain an expressionfor the correlators of fluctuations about the deterministic fixed-point of the system. This procedure is similar to thatused in [57–59]. These correlators are the basis for the theory results presented in the main text. Finally, we detailhow this analysis can be used in conjunction with the model detailed in Section S1 to produce the results in the maintext. A. Path-integral approach
We begin our analysis by defining a stochastic process in terms of the scaled variables q αj,t = n αj,t /N (we use a morecompact notation for the numbers of molecules n αj,t ≡ n αj ( t ) than in the main text). Here, N (defined in SectionS1) characterises the typical number of particles per cell and is sometimes referred to as the system size [70]. Tosimplify matters, we discretise time in steps of length ∆. The continuum limit is later recovered by taking ∆ → j , particles may be created or annihilated immediately and/or at one other future time t + τ . We define the number of reactions of type r which occur in the time interval t to t + ∆ and have an associateddelayed effect at t + τ by k j,r,τ,t . The number of particles of type α [where here α ∈ { ( p ) d, ( p ) h, ( m ) d, ( m ) h } ] whichare immediately created/annihilated in such a reaction is denoted ν αr and the number which are created/annihilatedat the later time t + τ is denoted w ατ,r .The stochastic process can therefore be written as q αj,t +∆ − q αj,t = (cid:88) r,τ ν αr k j,r,τ,t N + (cid:88) r,τ w ατ,r k j,r,τ,t − τ N . (S3)We wish to approximate this process with a set of stochastic differential equations. An elegant way to obtain thisapproximation is to start from a path-integral representation for the process, and then to perform an expansion ininverse powers of the system size (similar to that used by van Kampen [70]) within this representation. Using a path-integral approach, as opposed to a master equation, avoids the complications which arise due to the non-Markoviannature of the dynamics.We can write an expression for the probability of observing the system in a configuration n (which represents setof particle numbers of all types and locations) at time t as follows P ( n , t ) = 1 N (cid:90) (cid:89) j,α t (cid:89) t (cid:48) =0 dq αj,t (cid:48) δ (cid:20) N n − q ( t ) (cid:21) P (cid:0) { q αj,t (cid:48) } (cid:1) , (S4)where P (cid:0) { q αj,t (cid:48) } (cid:1) is the probability density of observing a particular trajectory (realisation of the system) { q αj,t (cid:48) } and δ ( · ) is the Dirac delta-function. Eq. (S4) is merely a statement that the probability of the system being in state n attime t is equal to the sum of the probabilities of observing any one of a set of paths which lead to the system beingin state n at time t .Constraining the trajectory to obey Eq. (S3) and then rewriting the Dirac delta functions as integrals of complexexponentials, one obtains P ( n , t ) = 1 N (cid:88) { k } (cid:90) (cid:89) j,α t (cid:89) t (cid:48) =0 dq αj,t (cid:48) δ (cid:20) N n − q ( t ) (cid:21) (cid:89) j,α,τ,t (cid:48) δ (cid:32) q αj,t (cid:48) +∆ − q αj,t (cid:48) − (cid:88) r,τ ν αr k j,r,τ,t (cid:48) N − (cid:88) r,τ w ατ,r k j,r,τ,t (cid:48) − τ N (cid:33) P ( { k } )3= 1 N (cid:88) { k } (cid:90) (cid:89) j,α t (cid:89) t (cid:48) =0 dq αj,t (cid:48) dp αj,t (cid:48) π δ (cid:20) N n − q ( t ) (cid:21) e i (cid:80) j,α,t (cid:48) p αj,t (cid:48) (cid:18) q αj,t (cid:48) +∆ − q αj,t (cid:48) − (cid:80) r,τ ναr kj,r,τ,t (cid:48) N − (cid:80) r,τ wατ,rkj,r,τ,t (cid:48)− τN (cid:19) P ( { k } ) , (S5)where P ( { k } ) is the probability of observing a particular set of creation/annihilation events { k } .We presume that each reaction event is independent such that P ( { k } ) = (cid:89) j,r,τ,t P ( k j,r,τ,t ) . (S6)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,τ,t be the probability per unit time squared (or ‘rate’) of a reaction of type r occurring at position j at atime between t and t + ∆ with associated delay between τ and τ + ∆. These rates may be dependent on the numbersof particles n . We presume that the numbers of reactions triggered in the interval t to t + ∆, k j,r,τ,t , are Poissonrandom variables with mean W j,r,τ,t ∆ . This involves the approximation that the reaction rates do not change overthe course of the small time step ∆. This is similar to the approximation made for the so-called tau-leaping variantof the Gillespie algorithm [71]. Making these assumptions, the distributions of the { k } are given by P ( k j,r,τ,t ) = (cid:0) ∆ W j,r,τ,t (cid:1) k j,r,τ,t k j,r,τ,t ! e − ∆ W j,r,τ,t . (S7)We note that in our system, the local reaction rate W j,r,τ,t can depend upon the number of particles in the adjacentcells. The sums over { k } in Eq. (S5) can be evaluated by observing that (cid:88) k j,r,τ,t (cid:0) ∆ W j,r,τ,t (cid:1) k j,r,τ,t k j,r,τ,t ! e − ∆ W j,r,τ,t e − i (cid:80) α p αj,t ναr kj,r,τ,tN − i (cid:80) α p αj,t + τ wατ,rkj,r,τ,tN = exp (cid:26) − ∆ W j,r,τ,t (cid:20) − e − i (cid:80) α ναr pαj,t + wατ,rpαj,t + τN (cid:21)(cid:27) . (S8)We finally arrive at the following expression for the probability distribution P ( n , t ) P ( n , t ) = 1 N (cid:90) (cid:89) j,α t (cid:89) t (cid:48) =0 dq αj,t (cid:48) dp αj,t (cid:48) π δ (cid:20) N n − q ( t ) (cid:21) e i (cid:80) j,α,t (cid:48) p αj,t (cid:48) ( q αj,t (cid:48) +∆ − q αj,t (cid:48) ) × exp − (cid:88) j,r,τ,t (cid:48) ∆ W j,r,τ,t (cid:48) (cid:20) − e − i (cid:80) α ναr pαj,t (cid:48) + wατ,rpαj,t (cid:48) + τN (cid:21) . (S9)Expanding the exponentials in Eq. (S9) in powers of 1 /N and truncating the series at next-to-leading order, oneobtains the following expression, P ( n , t ) = 1 N (cid:90) (cid:89) j,α t (cid:89) t (cid:48) =0 dq αj,t (cid:48) dp αj,t (cid:48) π δ (cid:20) N n − q ( t ) (cid:21) e i (cid:80) j,α,t (cid:48) p αj,t (cid:48) ( q αj,t (cid:48) +∆ − q αj,t (cid:48) ) × exp (cid:40) − (cid:88) j,r,τ,t (cid:48) ∆ W j,r,τ,t (cid:48) (cid:34) iN (cid:88) α (cid:0) ν αr p αj,t (cid:48) + w ατ,r p αj,t (cid:48) + τ (cid:1) + 12 N (cid:88) α,α (cid:48) ,j (cid:48) ,t (cid:48) ,τ (cid:48) δ τ,τ (cid:48) δ j,j (cid:48) δ t (cid:48) ,t (cid:48)(cid:48) (cid:0) ν αr p αj,t (cid:48) + w ατ,r p αj,t (cid:48) + τ (cid:1) (cid:16) ν α (cid:48) r p α (cid:48) j (cid:48) ,t (cid:48)(cid:48) + w α (cid:48) τ,r p α (cid:48) j (cid:48) ,t (cid:48)(cid:48) + τ (cid:48) (cid:17) (cid:35)(cid:41) , (S10)where δ l,l (cid:48) is the Kronecker delta. Eq. (S10) is similar to the Martin-Siggia-Rose-Janssen-de Dominicis (MSRJD)functional integral [69]. This result can be compared with the path integral for an SDE of an appropriate form. Thecorresponding path integral expression for the general SDE dn αj,t dt = F αj,t ( n t ) + ξ αj,t , (cid:104) ξ αj,t ξ α (cid:48) j (cid:48) ,t (cid:48) (cid:105) = Σ αα (cid:48) j,j (cid:48) ,t,t (cid:48) , (S11)is given by [57–59] P ( n , t ) = 1 N (cid:90) (cid:89) j,α t (cid:89) t (cid:48) =0 dq αj,t (cid:48) dp αj,t (cid:48) π δ (cid:20) N n − q ( t ) (cid:21) e i (cid:80) j,α,t (cid:48) p αj,t (cid:48) ( q αj,t (cid:48) +∆ − q αj,t (cid:48) ) × exp (cid:34) − i ∆ N (cid:88) α,j,t (cid:48) p αj,t (cid:48) F αj,t (cid:48) − ∆ N (cid:88) α,α (cid:48) ,j,j (cid:48) ,t (cid:48) ,t (cid:48)(cid:48) p αj,t (cid:48) p α (cid:48) j (cid:48) ,t (cid:48)(cid:48) Σ αα (cid:48) j,j (cid:48) ,t (cid:48) ,t (cid:48)(cid:48) (cid:35) . (S12)From Eq. (S10) one can then read off the following effective SDEs, which are good approximations of the dynamicsfor large but finite N , by taking the limit ∆ → dn αj,t dt = (cid:88) r W j,r,t ν αr + (cid:88) r (cid:90) ∞ W j,r,τ,t − τ w αr,τ dτ + ξ αj,t , (S13)where W j,r,t = (cid:82) ∞ W j,r,τ,t dτ , and where the correlators of the stochastic noise variables are given by (cid:104) ξ αj,t ξ α (cid:48) j (cid:48) ,t (cid:48) (cid:105) = δ ( t − t (cid:48) ) δ j,j (cid:48) (cid:34)(cid:88) r ν αr ν α (cid:48) r W j,r,t + (cid:88) r (cid:90) ∞ w αr,τ w α (cid:48) r,τ W j,r,τ,t − τ dτ (cid:35) + δ j,j (cid:48) (cid:34)(cid:88) r W j,r,t − t (cid:48) ,t w αr,t − t (cid:48) ν α (cid:48) r + (cid:88) r W j,r,t (cid:48) − t,t w α (cid:48) r,t (cid:48) − t ν αr (cid:35) . (S14) B. 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,τ,t = W j,r,t K r ( τ ). That is, the delay time τ is drawn from a distribution K r ( τ ) which is independent of j , t and n αj,t . Eq. (S13) then becomes dn αj,t dt = (cid:88) r W j,r,t ν αr + (cid:88) r (cid:90) ∞ W j,r,t − τ K r ( τ ) w αr,τ dτ + ξ αj,t . (S15)If we consider small fluctuations about the fixed point of the deterministic system δ αj,t = n αj,t − ¯ n α [as in Eq. (1)],Eq. (S15) may be approximated by dδ αj,t dt = (cid:88) j (cid:48) ,α (cid:48) J α,α (cid:48) j,j (cid:48) δ α (cid:48) j (cid:48) ,t + (cid:88) j (cid:48) ,α (cid:48) (cid:90) ∞−∞ L α,α (cid:48) j,j (cid:48) ,τ δ α (cid:48) j (cid:48) ,t − τ dτ + ξ αj,t , (S16)where J α,α (cid:48) j,j (cid:48) = (cid:32)(cid:88) r ∂W j,r,t ∂n α (cid:48) j (cid:48) ,t ν αr (cid:33) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( n j (cid:48) ,t =¯ n ) , (S17)and L α,α (cid:48) j,j (cid:48) ,τ = (cid:32)(cid:88) r ∂W j,r,t − τ ∂n α (cid:48) j (cid:48) ,t − τ K r ( τ ) θ ( τ ) w αr,τ (cid:33) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( n j (cid:48) ,t − τ =¯ n ) , (S18)where θ ( τ ) is the Heaviside function. We note that for the systems studied in the main text, J α,α (cid:48) j,j (cid:48) and L α,α (cid:48) j,j (cid:48) ,τ arenon-zero for j (cid:54) = j (cid:48) 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 theevaluation of the correlators Eq. (S14). That is, we evaluate (cid:104) ξ αj,t ξ α (cid:48) j (cid:48) ,t (cid:48) (cid:105) at n αj,t = ¯ n α . As result, what was multiplicativenoise is now treated as additive, simplifying the calculation.5Carrying out a temporal Fourier transform in Eq. (S16), one obtains (cid:88) α (cid:48) ,j (cid:48) (cid:104) iωδ α,α (cid:48) δ j,j (cid:48) − J α,α (cid:48) j,j (cid:48) − ˆ L α,α (cid:48) j,j (cid:48) ,ω (cid:105) ˆ δ α (cid:48) j (cid:48) ,ω = ˆ ξ αj,ω . (S19)Eq. (S19) can be rewritten in matrix form as M − ω ˆ δ ω = ˆ ξ ω , (S20)where the different elements of the vector correspond to the different types of particle at the different cell sites. Writingthe correlation matrix of the noise variables as Σ ω = (cid:104) ˆ ξ ω ˆ ξ † ω (cid:105) one finally obtains the following result for the matrix ofthe correlators of the fluctuations C ( ω ) ≡ (cid:104) ˆ δ ω ˆ δ † ω (cid:105) = M ω Σ ω M † ω . (S21)The diagonal elements of this matrix correspond to the power spectrum of fluctuations in Eq. (2), i.e. P αj ( ω ) = C α,αj,j ( ω ). The off-diagonal elements allow one to calculate the phase lag as stated in Eq. (3), that is C α,α (cid:48) j,j (cid:48) ( ω ) = (cid:104) ˆ δ αj,ω ˆ δ α (cid:48) (cid:63)j (cid:48) ,ω (cid:105) .When the problem is translationally invariant, as is the case when the system parameters are the same in all cells, J α,α (cid:48) j,j (cid:48) and L α,α (cid:48) j,j (cid:48) ,τ are functions of j − j (cid:48) only. One can then further simplify matters by carrying out a Fourier transformwith respect to position as well. One obtains (cid:88) α (cid:48) (cid:104) iωδ α,α (cid:48) − ˜ J α,α (cid:48) k − ˆ˜ L α,α (cid:48) k,ω (cid:105) ˆ˜ δ α (cid:48) k,ω = ˆ˜ ξ αk,ω , (S22)and from this M − k,ω ˆ˜ δ k,ω = ˆ˜ ξ k,ω , (S23)where the different elements of the vector now correspond only to the different species. We then arrive at a similarexpression to Eq. (S21), but where the matrix dimension is reduced by a factor of L (the number of cells) C ( k, ω ) ≡ (cid:104) ˆ˜ δ k,ω ˆ˜ δ † k,ω (cid:105) = M k,ω Σ k,ω M † ω . (S24)Here, the diagonal elements correspond to the power spectra in Eq. (6) and the off-diagonal elements allow one tocalculate the phase lag using Eq. (3) of the main text. C. Application to the model gene regulatory network
Using Eqs. (S13) and (S14), the individual-based system given by Eq. (S1) can be approximated by stochasticdifferential equations of the form dn (p)h j,t dt = a h n (m)h j,t − τ (p)h − b h n (p)h j,t + ξ (p)h j,t ,dn (p)d j,t dt = a d n (m)d j,t − τ (p)d − b d n (p)d j,t + ξ (p)d j,t ,dn (m)h j,t dt = f h j (cid:0) n t − τ (m)h (cid:1) − c h n (m)h j,t + ξ (m)h j,t ,dn (m)d j,t dt = f d j (cid:0) n t − τ (m)d (cid:1) − c d n (m)d j,t + ξ (m)d j,t . (S25)The quantities involved in these equations are given in Section S1, with the exception of the Gaussian noise terms ξ αj,t ,which encapsulate the intrinsic noise due to the stochastic and individual-based nature of the gene regulatory system.The deterministic trajectories ¯ n αj,t , examples of which are given in Figs. 2 and 3, are found by setting the noise terms ξ αj,t in Eq. (S25) to zero and numerically integrating the resulting ordinary differential equations. Setting the noise6terms to zero effectively approximates the system as being infinitely large, so that any fluctuations are negligible incomparison to the mean particle numbers. This is not appropriate in cases in which intrinsic noise significantly affectsthe dynamics (as is the case for the examples we look at). This is exemplified by the stark disagreement between thedeterministic and individual-based simulations in Figs. 2 and 3.When the deterministic (infinite) system has reached a fixed point, the noise terms ξ αj,t in Eq. (S25) can be takento have the following correlators within the linear-noise approximation, (cid:104) ξ (p)h j,t ξ (p)h j (cid:48) ,t (cid:48) (cid:105) = δ ( t − t (cid:48) ) δ j,j (cid:48) (cid:104) a h ¯ n (m)h + b h ¯ n (p)h (cid:105) , (cid:104) ξ (p)d j,t ξ (p)d j (cid:48) ,t (cid:48) (cid:105) = δ ( t − t (cid:48) ) δ j,j (cid:48) (cid:104) a d ¯ n (m)d + b d ¯ n (p)d (cid:105) , (cid:104) ξ (m)h j,t ξ (m)h j (cid:48) ,t (cid:48) (cid:105) = δ ( t − t (cid:48) ) δ j,j (cid:48) (cid:104) f h ( ¯n ) + c h ¯ n (m)h (cid:105) , (cid:104) ξ (m)d j,t ξ (m)d j (cid:48) ,t (cid:48) (cid:105) = δ ( t − t (cid:48) ) δ j,j (cid:48) (cid:104) f d ( ¯n ) + c d ¯ n (m)d (cid:105) , (S26)where barred quantities are evaluated at the deterministic fixed point (i.e. the solution to Eqs. (S25) with the noiseterms and the time derivatives set to zero). All inter-species cross-correlators are zero, due to the fact that no onereaction 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 M ω or M k,ω , defined in Eqs. (S20) or (S23) respectively, by performing the linear-noise approximation on Eq. (S25), asdetailed in Section S2 B. Using the correlators in Eq. (S26), one is then able to evaluate the matrix of correlators C ( ω ) or C ( k, ω ), defined in Eqs. (S21) or (S24) respectively, which contains all the information one needs to find thepower spectra of fluctuations and/or the phase lag (as is also described in Section S2 B).As an example, we provide expressions that allow one to calculate the matrix M k,ω in the case that the system isspatially homogeneous. In Eq. (S22), the matrix ˜ J αα (cid:48) k is given by˜ J k = − b h − b d − c h
00 0 0 − c d , (S27)and the matrix ˆ˜ L αα (cid:48) k,ω is given byˆ˜ L k,ω = a h e iωτ (p)h
00 0 0 a d e iωτ (p)d (cid:16) f h h e iωτ (m)h (cid:17) (cid:16) f h d e iωτ (m)h ( d (+) e − ik + d ( − ) e ik ) (cid:17) (cid:16) f d h e iωτ (m)d (cid:17) (cid:16) f d d e iωτ (m)d ( d (+) e − ik + d ( − ) e ik ) (cid:17) , (S28)where the entries of the matrix are in the order α = { (p)h , (p)d , (m)h , (m)d } and where we have f αh = − k α N n ( p ) h ( n ( p ) h + ¯ n ( p ) h ) (cid:34) r α h + r α hd ¯ n ( p ) d n ( p ) d + ¯ n ( p ) d (cid:35) ,f αd = k α N n ( p ) d ( n ( p ) d + ¯ n ( p ) d ) (cid:34) r α d + r α hd n ( p ) h n ( p ) h + ¯ n ( p ) h (cid:35) . (S29) S3. LINEAR STABILITY ANALYSIS
We have derived expressions for the power spectrum of fluctuations in the individual-based system about thedeterministic fixed point. Such an analysis presumes the stability of this fixed point. To determine whether or notsuch a stable fixed point exists for a particular parameter set, we perform a linear stability analysis of the deterministicsystem.We begin with the linearised expression for the time-evolution of small deviations about the deterministic fixedpoint Eq. (S16), but with the noise term removed dδ αj,t dt = (cid:88) j (cid:48) ,α (cid:48) J α,α (cid:48) j,j (cid:48) δ α (cid:48) j (cid:48) ,t + (cid:88) j (cid:48) ,α (cid:48) (cid:90) ∞ L α,α (cid:48) j,j (cid:48) ,τ δ α (cid:48) j (cid:48) ,t − τ dτ. (S30)7For simplicity, we now assume that the delay kernels are all Dirac-delta functions [i.e. K r ( τ ) = δ ( τ − τ r )] such that (cid:88) j (cid:48) ,α (cid:48) (cid:90) ∞ L α,α (cid:48) j,j (cid:48) ,τ δ α (cid:48) j (cid:48) ,t − τ dτ = (cid:88) j (cid:48) ,α (cid:48) ,r H α,α (cid:48) j,j (cid:48) ,r δ α (cid:48) j (cid:48) ,t − τ r , (S31)where H α,α (cid:48) j,j (cid:48) ,r = (cid:32) ∂W j,r,t − τ ∂n α (cid:48) j (cid:48) ,t − τ w αr (cid:33) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( n j (cid:48) ,t − τ =¯ n ) . (S32)We have further supposed that the number of particles produced by a delayed reaction w αr is independent of the delaytime τ r .We then suppose that the deviation away from the fixed point evolves in an exponential way. That is we proposethe ansatz δ αj,t = δ α (0) j e λt . (S33)Upon substitution into Eq. (S30), one obtains (cid:88) j (cid:48) ,α (cid:48) (cid:34) λδ j,j (cid:48) δ α,α (cid:48) − J α,α (cid:48) j,j (cid:48) − (cid:88) r H α,α (cid:48) j,j (cid:48) ,r e λτ r (cid:35) δ α (cid:48) (0) j (cid:48) = 0 . (S34)We note that if we had not assumed a delta function for the delay kernel K r ( τ ), the exponential factors in Eq. (S35)would instead be replaced by L τ { K r ( τ ) } ( − λ ), the Laplace transform of the delay kernel evaluated at − λ . Eq. (S34)can be rewritten more succinctly as the following matrix equation (cid:34) λ − J − (cid:88) r H r e λτ r (cid:35) δ (0) = 0 . (S35)In order for the solutions to this equation to be non-trivial, the determinant of the object multiplying δ (0) in Eq. (S35)must be equal to zero. This gives rise to an ‘eigenvalue’ equation for λ . Due to the exponential terms involving λ which arise from the delays, this equation is not analytically tractable and must be solved numerically.Typically, one finds many possible solutions λ to the eigenvalue equation for any one set of system parameters, andthus the full solution to Eq. (S30) is a linear combination of solutions of the form in Eq. (S33). If any one of theeigenvalues λ has a positive real part the fixed point is unstable. Else, we say that it is (linearly) stable. In order toproduce the dotted lines in Figs. 5 and 6, one finds the eigenvalues λ as a function of the system parameters. Theinstability line divides regions in parameter space for which all the eigenvalues have negative real parts from thoseregions for which some eigenvalues have positive real parts. The imaginary part of the least stable eigenvalue for eachparameter set is plotted as a green line in Figs. 5a and 6a.Once one has found the eigenvalues λ , one can then also go on to solve (numerically) for the eigenvectors δ (0)(0)