A unified mechanism for spatiotemporal patterns in somitogenesis
AA unified mechanism for spatiotemporal patterns in somitogenesis
Chandrashekar Kuyyamudi,
1, 2
Shakti N. Menon, and Sitabhra Sinha
1, 2 The Institute of Mathematical Sciences, CIT Campus, Taramani, Chennai 600113, India Homi Bhabha National Institute, Anushaktinagar, Mumbai 400 094, India (Dated: March 30, 2020)Somitogenesis, the process of body segmentation during embryonic development, exhibits a keyset of features that is conserved across all vertebrate species despite differences in the detailedmechanisms. Prior to the formation of somites along the pre-somitic mesoderm (PSM), periodicexpression of clock genes is observed in its constituent cells. As the PSM expands through theaddition of new cells at its posterior, the oscillations in the cells closer to the anterior cease andeventually lead to formation of rostral and caudal halves of the somites. This pattern formation isbelieved to be coordinated by interactions between neighboring cells via receptor-ligand coupling.However, the mechanism underlying the transition from synchronized oscillations to traveling wavesand subsequent arrest of activity, followed by the appearance of polarized somites, has not yet beenestablished. In this paper we have proposed a unified mechanism that reproduces the sequence ofdynamical transitions observed during somitogenesis by combining the local interactions mediatedvia Notch-Delta intercellular coupling with global spatial heterogeneity introduced through a mor-phogen gradient that is known to occur along the anteroposterior axis of the growing PSM. Ourmodel provides a framework that integrates a boundary-organized pattern formation mechanism,which uses positional information provided by a morphogen gradient, with the coupling-mediatedself-organized emergence of collective dynamics, to explain the processes that lead to segmentation.
I. INTRODUCTION
Somitogenesis is the process of formation of somites,which are the modular building blocks of all vertebratebodies [1–3]. Somites compose bilaterally symmetric seg-ments that are formed in the paraxial, or pre-somitic,mesoderm (PSM) of developing embryos as the bodyaxis itself elongates [4]. Analogous processes have beenimplicated in the body segmentation of some inverte-brates [5, 6]. Although there is great variability acrossspecies in terms of the number of somites, the meansize of a somite and the duration over which they areformed, nonetheless a conserved set of features character-izing somitogenesis is seen across these species [7]. A gen-eral conceptual model for explaining these core featuresis provided by the Clock and Wavefront (CW) frameworkproposed by Cooke and Zeeman in 1976 [8], wherein thePSM was assumed to comprise cellular oscillators (clocks)which are each arrested at their instantaneous state ofactivity upon encountering a wavefront that moves fromthe anterior to posterior of the PSM [9–16]. This pro-vides a route for translating a sequence of temporal ac-tivity into spatial patterns [17]. In order to constructan explicit mechanism embodying the CW framework,we need to disaggregate its components that operate atdifferent scales, namely, (i) the cellular scale at whichoscillations occur, (ii) the inter-cellular scale at whichcontact-mediated signaling takes place, and (iii) the scaleof the PSM across which morphogen gradients form andact as the environment that could modulate the inter-cellular interactions. This resonates with the proposalof Oates [18] to view the CW framework as a three-tierprocess. In the bottom tier, we observe oscillations atthe level of a single cell in the PSM, arising from theperiodic expression of clock genes [19–26]. The middle tier describes the mechanism by which the cellular oscil-lators coordinate their activity with that of their neigh-bors. This occurs through juxtacrine signaling broughtabout by interactions between Notch receptors and Deltaligands [22, 27–38]. Indeed, several earlier models haveexplored the role of Notch-Delta coupling in bringingabout robust synchronization between the oscillators [39–43]. Finally, processes that bring about the slowing down(and eventual termination) of the oscillations [19, 32, 44],and the subsequent differentiation of the cells into ros-tral and caudal halves of the somites, constitute the toptier. The model proposed here integrates these differ-ent scales by investigating genetic oscillators interactingvia Notch-Delta coupling whose strength is modulated ina position-dependent manner due to a morphogen con-centration gradient along the anteroposterior (AP) axisof the PSM. This provides an unified framework for ex-plaining the dynamical transitions observed during somi-togenesis.As the PSM expands along the AP axis through theaddition of new cells at the tail [4], it is reasonable torestrict our attention to the process of somitogenesistaking place along a one-dimensional array of coupledcells in the PSM. From the perspective of our model-ing which explicitly investigates the role of morphogengradient in coordinating somite formation, the array ofcells is considered to be aligned along the AP axis, as thevarious morphogens that are expressed along the PSMare known to form concentration gradients along thisaxis [45]. These primarily include molecules belonging tothe FGF [46, 47], RA [48, 49] and Wnt families [50–52].Even though the role of gradients on the overall dynamicshas been explored [53–56], there is to date no consensusas to the explicit mechanism through which they con-tribute to somite formation. Our model demonstrates a r X i v : . [ q - b i o . T O ] M a r FIG. 1.
The key dynamical features of somitogenesis(top) that are reproduced in our mathematical model(bottom). (a) Schematic diagram depicting the zebrafishpre-somitic mesoderm (PSM) and the dynamical states ex-hibited by the cells over the course of somitogenesis. Thebox outlined with broken lines represents the posterior regionwhich exhibits a spatial gradient in morphogen concentra-tion, represented by the exponentially decaying profile shownabove the box. The dynamics of the cells at the tail (regionIV) are characterized by exact synchronization of the periodicvariation in gene expression. As the cells move towards theanterior, they first exhibit travelling waves of gene expression(III), followed by arrest of the oscillations (II). Eventually thecells differentiate (I) into alternating bands corresponding torostral and caudal halves of the mature somites. The figures inthe insets at the right display typical time series in each of theregimes I-IV of the inhibitor gene expression for two neigh-boring cells coupled to each other. (b) Schematic diagramdescribing the interaction between two cells via Notch-Deltacoupling. In general, each cell has Notch receptors, as well as,Delta ligands that bind to them. The cis form of the bindingleads to the loss of receptors and ligands without resulting inany downstream signaling, whereas trans
Notch-Delta bind-ing gives rise to cleavage of the Notch Intra-Cellular Domain(NICD). The latter acts as a transcription factor (TF) for thedownstream activator ( x ) and inhibitor ( y ) genes which arethe essential constituents of each cellular oscillator. that if the morphogen gradient is considered to regulatethe strength with which adjacent cells interact, it canlead to qualitatively different kinds of dynamics alongthe PSM in a threshold-dependent manner [Fig. 1 (a)].Unlike the conventional boundary-organized pattern for-mation paradigm [57], here the morphogen gradient doesnot determine the cell fate so much as affect the local dy-namics such that, at any point in time, there can be qual- itatively different types of dynamics occurring in separatelocations in the PSM. Thus, our results help address anopen question as to how morphogen gradients influencethe collective dynamics of the cellular oscillators, leadingto their eventual fates as rostral and caudal halves of themature polarized somites. II. METHODSA. Modeling genetic oscillators interacting viaNotch-Delta coupling
Several experiments have established that the cellsin the presomitic mesoderm (PSM) have “clock” geneswhose expression levels oscillate [19–24, 26]. In ourmodel, we consider a generic two component genetic os-cillator comprising an activator gene ( x ), which upregu-lates its own expression, as well as that of an inhibitorgene ( y ), which suppresses the expression of the activatorgene [58]. The dynamics of this two-component oscilla-tor can be expressed in terms of a pair of rate equationsdescribing the change in concentrations of the proteinproducts X and Y of genes x and y , respectively. Themodel parameters are chosen such that X and Y exhibitlimit cycle oscillations.As the communication between cells in the PSM is cru-cial in mediating their collective behavior during somito-genesis, we couple the dynamics of the clock genes ofneighboring cells. Experiments have established the roleof the Notch-Delta juxtacrine signaling pathway in me-diating the interaction between cells that are in physicalcontact with each other [22, 27–36]. In general, each cellhas both Notch receptors as well as Delta ligands on theirsurface. A Notch receptor on cell i which is bound to aDelta ligand belonging to a neighboring cell j (i.e., trans binding) leads to the cleavage of the Notch intracellulardomain (NICD) that will act as transcription factor fordownstream genes in cell i [59, 60]. In our model we as-sume that the NICD upregulates the expression of boththe clock genes, while the gene products X, Y suppressthe production of Delta ligands by the cell [39, 41]. Wedescribe the Notch-Delta signaling mechanism throughthe coupled dynamics of (i) the free (unbound) Notch re-ceptor concentration ( N ), (ii) the free Delta ligand ( D )and (iii) the NICD which is released as a result of trans binding ( N b ). Thus, the dynamics of a cell i coupled toits neighbors through Notch-Delta signaling is describedby the following set of equations: dX i dt = a + bX i + f N bi X i + Y i + N bi − cX i , (1) dY i dt = eX i + g i N bi X i + N bi − Y i , (2) dN i dt = β N − γN i − D i N i − D tr N i k , (3) dN bi dt = D tr N i k − µN bi , (4) dD i dt = β D h ( X i + Y i ) − γD i − D i N i − D i N tr k . (5)Here the terms D tr and N tr are the mean values of D j and N j over all neighboring cells j to which i is coupledthrough trans -binding. While the values of the modelparameters can, in general, vary across cells, we restrictour attention to the variation of the coupling parame-ter g (subscripted with the cell index in Eqns. (1-5)),which determines the strength of upregulation of y by theNICD ( N b ). This allows us to investigate the role of spa-tial heterogeneity imposed by the gradient of morphogenconcentration along the anteroposterior (AP) axis of thePSM. B. Dynamical evolution of the morphogen gradient
Our model focuses on the behavior of a contiguoussegment of cells of length (cid:96) in the PSM with a mor-phogen source located at its anterior. If the strengthof the source is constant in time, it would have resultedin the gradient becoming progressively less steep as thePSM expands. This dilution is countered by the net in-crease in the strength of the source through the secretionof morphogen by the newly matured somites. Thus, asthe PSM expands due to addition of cells at the poste-rior tail, we can view the segment under considerationas effectively flowing up a morphogen gradient along theAP axis [18]. We choose a segment of N cells with aspatial extent (cid:96) that is initially located ( t = 0) at theposterior end of the PSM, i.e., at the lower end of thegradient, and follow its evolution upto a time period T .As mentioned above, the effect of the varying morphogenconcentration on the dynamics of the cells is introducedvia the coupling parameter g . Specifically, we assumethat the value of g at each site is proportional to thecorresponding morphogen concentration, yielding an ex-ponentially decaying gradient of g across the AP axis: g i ( t ) = g min exp ( λ g x i ( t )). The steepness of the gradi-ent is quantified by λ g , which is a function of T , as wellas g max and g min , which are the values of g at the an-terior and posterior ends of the PSM, respectively, viz., λ g = ln( g max /g min ) /T . We assume that the effective flowof the segment of cells along the AP axis occurs at an uni-form rate. This can be taken to be unity without loss ofgenerality by appropriate choice of time unit. Thus, theinstantaneous position x i ( t ) along the gradient of the i thcell in the segment is given by x i ( t ) = t + ( (cid:96)/N )( i − x ( t = 0) = 0. III. RESULTS
In our simulations, we have considered the PSM tocomprise cells, each of which exhibits oscillating gene ex- pression. We assume a minimal model for the geneticoscillator consisting of two clock genes, one activatoryand the other inhibitory, whose products correspond tofate determining proteins (see Methods). The oscillationsof neighboring cells influence each other through Notch-Delta inter-cellular coupling [Fig. 1 (b)]. We have explic-itly verified that incorporating delay in the contact medi-ated signaling does not alter our results qualitatively (seeSI). Considering the simplest setting, viz., a pair of adja-cent cells, which allows us to investigate the effect of cou-pling on the collective dynamics, we see from Fig. 2 thatthe system can exhibit a wide range of spatio-temporalpatterns. These can be classified systematically throughthe use of quantitative measures (see Supplementary In-formation for details). We focus on the patterns thatcan be immediately interpreted in the context of somi-togenesis: (i) Inhomogeneous Steady States (ISS), (ii)Homogeneous Steady States (HSS), (iii) Anti-Phase Syn-chronization (APS) and (iv) Exact Synchronization (ES)[shown schematically as insets of Fig. 1 (a)]. The range ofvalues of the coupling-related parameters f , g and h overwhich these patterns are observed in a pair of coupledcells are shown in Fig. 2. The parameters f and g governthe strength with which the Notch intra-cellular domain( N b ) regulates the activatory and inhibitory clock genes,respectively, while h is related to the intensity of repres-sion of the Delta ligand ( D ) by each of the clock genes.As inter-cellular coupling is believed to be responsiblefor the synchronized activity of cells in the initial stageof somitogenesis [61], we note that the dynamical regimecorresponding to ES occurs for low g and intermediatevalues of f , with the region increasing in size for larger h . Spatial variation in these coupling parameters acrossthe PSM can arise through heterogeneity in the underly-ing morphogen concentrations. It is known that the mor-phogens RA, Wnt and FGF are differentially expressedalong the PSM, exhibiting monotonically varying con-centration gradients having peaks at the posterior (forFGF and Wnt) or anterior (for RA) ends [51, 62, 63].Experiments on several vertebrate species have shownthat high concentrations of RA initiate differentiation,while increased levels of Wnt and FGF, which are knownto play a role in sustaining oscillations in the expres-sion of the clock genes, impede the formation of ma-ture somites [46, 49, 64, 65]. Thus, as our goal is toexplicate the mechanisms driving the differentiation ofindividual cells into anterior and posterior halves of themature somites, we focus on the role of RA on the col-lective dynamics of cells in the PSM. We incorporate theeffect of this morphogen in the spatial variation of thecoupling parameter g which is assumed to exponentiallydecay from the anterior to the posterior end of the do-main. Such a profile will naturally arise if the morphogendiffuses from a source located at the anterior and is de-graded at a constant rate across space [66–68]. We notethat qualitatively similar results are obtained if gradientsare introduced in either f or h instead of g . The effect of FIG. 2.
Transitions between patterns representing dif-ferent stages in somitogenesis seen in the collectivedynamics of a pair of cells on varying the parametersgoverning the strength of Notch-Delta coupling inour model. (a) Schematic diagram of the three-dimensionalspace spanned by the coupling parameters ( f, g, h ), scaled bythe relevant kinetic parameters of the individual oscillators.Note the logarithmic scale used for the ranges of the param-eters. (b) The variation with g of the relative frequency ofoccurrence of patterns belonging to each of six distinct cat-egories, viz., Exact Synchronization (ES), Anti-Phase Syn-chronization (APS), Homogeneous Steady States (HSS), In-homogeneous Steady States (ISS), Inhomogeneous In-phaseSynchronization (IIS) and Others. The values of the param-eters f/b and h are fixed at 0 .
25 and 4, respectively. (c-n)The most commonly occurring dynamical patterns (i.e., ob-tained for >
50% of all initial conditions used) that are seenfor different values of f/b and h (varying over four order ofmagnitude) for 12 equally spaced values of log( g/e ) between − .
89 (panel c) and 0 .
36 (panel n). While HSS is the mostcommon pattern seen over this range of parameter values,focusing on how the occurrence frequency of ES varies with f, g, h indicates where a transition from synchronization totime-invariant behavior may be achieved.
Wnt and FGF are implicitly accounted for by the endo-geneously oscillating system of two coupled clock genesthat we consider in our model.Introducing heterogeneity through the coupling pa-rameter in the pair of adjacent cells considered earlier, a |g −g |( × −2 ) < g > b −g |( × −2 ) < g > c ∆ d ∆ e A τ P t (a.u.)Y 60 1203.03.6 t (a.u.)A f
60 1202.12.6 t (a.u.) τ P g FIG. 3.
Incorporating a morphogen gradient by vary-ing the Notch-Delta coupling parameter g of adja-cent oscillators reproduces the temporal sequence ofpatterns observed during somitogenesis. (a) Schematicrepresentation of the spatial variation of g , resulting fromthe different concentrations of a morphogen sensed by neigh-boring cells. The steepness of the gradient is quantified by | g − g | , the difference in the values g , for the two oscil-lators, while their location on the gradient is determined bythe mean (cid:104) g (cid:105) . (b) The most commonly occurring dynami-cal patterns (i.e., obtained for >
50% of all initial conditionsused) on varying | g − g | and (cid:104) g (cid:105) . The steady state (SS) re-gion, representing both HSS and ISS patterns, spans a largearea across parameter space. In this region, the states of theadjacent cells, which are characterized by the protein concen-tration Y , converge to fixed points Y , . (c) The differencebetween the steady state concentrations of the adjacent cells,∆ = | Y − Y | , increases with (cid:104) g (cid:105) as one effectively moves upthe morphogen gradient, while being relatively unaffected bythe steepness | g − g | . (d-g) The dynamical consequences of amorphogen gradient with an exponentially decaying concen-tration profile. (d) As a pair of coupled cells gradually moveupstream of the gradient, resulting in an increase of g , overtime t , their collective behavior (represented by Y ) convergesfrom initially synchronized oscillations (ES) to an inhomoge-neous steady state (ISS, characterized by finite values of ∆)at long times [see SI, Fig. S5]. The changes occurring in thesystem during the transition from oscillations to steady statebehavior (shaded region) can be quantitatively investigatedby focusing on how the amplitude A and period τ p of the os-cillations (see panel e) change over time. (f-g) As cells moveupstream of the morphogen gradient (corresponding to a pro-gression from posterior to anterior regions in the PSM), themodel exhibits decreasing A (f) and increasing τ p (g). This isconsistent with a key experimental observation, viz., shorten-ing and slowing of oscillations as cells approach the anteriorend of the PSM, during somitogenesis [18]. we observe that qualitatively similar spatio-temporal pat-terns to those observed in Fig. 2 are obtained (see Sup-plementary Information, Fig. S3). On varying the meanvalue of the coupling (cid:104) g (cid:105) = ( g + g ) /
2, which effectivelyrepresents the location of these cells on the PSM, andthe steepness of the gradient | g − g | [Fig. 3 (a)], therange of (cid:104) g (cid:105) over which ES is seen (corresponding to theregion proximal to the posterior end of the PSM) doesnot appear to change appreciably on increasing | g − g | [Fig. 3 (b)]. Above a critical value of (cid:104) g (cid:105) which is indepen-dent of the gradient, the activities of the cells are arrestedat Y and Y , respectively, with the gap ∆ = | Y − Y | becoming larger as we move towards the anterior end,corresponding to increasing (cid:104) g (cid:105) [Fig. 3 (c)].As explained in the Methods, over the course of devel-opment, the PSM expands through new cells being addedto its posterior end, such that the existing cells progres-sively encounter increasing values of the morphogen con-centration. Modeling this time-evolution as an effectiveflow of the segment of adjacent cells along the gradient in g , we observe that a transient phase of ES is followed bydesynchronization and subsequent attenuation of the os-cillations, eventually leading to a separation of the steadystates of the two cells [Fig. 3 (d)]. The gap ∆ betweenthe steady states increases with time, giving rise to a pro-nounced ISS state. Immediately preceding the arrest ofperiodic activity, we observe that the amplitude A andperiod τ P [Fig. 3 (e)] of the oscillations show charac-teristic changes that are in agreement with experimen-tal observations of somitogenesis. Specifically, the am-plitude decreases [Fig. 3 (f)] while the period increases[Fig. 3 (g)] with time, similar to what has been reportedin Refs. [18, 69].Having seen that a pair of contiguous cells can in-deed converge to markedly different steady state valuesof their clock gene expressions, we now investigate thegeneralization to a spatially extended segment of length (cid:96) in the growing PSM subject to a morphogen gradi-ent [varying from g min to g max as shown schematicallyabove Fig. 4 (a)]. As the principal variation of the mor-phogen concentration occurs along the AP axis of thePSM, we restrict our focus to a one-dimensional array ofcells aligned along this axis. As new cells are added tothe posterior of the PSM over time, the relative positionof the segment of cells under consideration shifts alongthe morphogen gradient from the posterior to the ante-rior. We note that had there been a temporally invariantsource of morphogen at the anterior, the expansion ofthe PSM would have resulted in a dilution of the gra-dient. However, newly matured somites at the anteriorend serve as additional sources of the morphogen overtime [70–72], thereby ensuring that each cell experiencesan exponential increase in the morphogen concentration(see SI). This is reproduced in our model by the seg-ment effectively flowing up the morphogen gradient (asdiscussed in Methods).As shown in Fig. 4 (a), for a range of values of the pa-rameters g max , g min and (cid:96) , the cells display a short-livedES pattern, which is followed by the development of aphase lag between adjacent cells [as can be seen in theinset of Fig. 4 (b)]. This is analogous to the appearance c e ll i nde x a
120 Y0 2 4 6 8 10 120 300510Y t (a.u.) b FIG. 4.
Collective dynamics of a cellular array re-sponding to an exponential gradient of morphogenconcentration reproduces the spatio-temporal evolu-tion of PSM activity seen during somitogenesis.
Thetransition from progenitor cells in the posterior (P) to matu-rity at the anterior (A) of a segment of length (cid:96) comprising N cells in the PSM viewed as a flow upstream (moving windowin the schematic on top) over a period of time T along theexponential profile of the parameter g , decaying from g max to g min , reflecting the morphogen gradient. Different cells inthe window sense different morphogen concentrations, whosevalues change over time. This is incorporated in terms of thetime-dependent gradient ∆ g ( t ), with g ( t ) and g N ( t ) beingthe values of the coupling parameter g at the anterior andposterior ends of the window at time t , respectively. The re-sulting change in the activity of a segment comprising N = 20coupled cells, as it moves from P to A, is shown in (a). Thesystem initially exhibits synchronized oscillations across thesegment but, as a consequence of the gradient, a phase lagdevelops between adjacent cells (as seen in the time seriesin panel b). This results in a wave-like propagation of thepeak expression from the anterior to the posterior end of thesegment in each cycle. Subsequently, the oscillations reducein intensity leading to a homogeneous steady state (HSS),a magnified view of the transition being shown in the insetof panel (b). Eventually, the system converges to an inho-mogeneous steady state (ISS), with adjacent cells attainingdifferent fates characterized by alternating high and low val-ues of the protein concentration Y [cells with odd and evenindices on the segment are shown using different colors inpanel (b)]. Note that the system exhibits the entire range ofpatterns shown in the schematic in Fig. 1 (a), with the tem-porally invariant spatial pattern observed at the anterior endof the PSM resembling that of alternating rostral and caudalhalves of somites. Results shown are for parameter values of g max = 15 . g min (= g ) = 1 . T = 300 . of a small phase difference between the pair of oscillatorsdescribed earlier, and manifests as a travelling wave thatpropagates along the segment. As the cells move furtherup the gradient, the oscillations subside and the segmentexhibits a transient homogeneous state, which is consis-tent with the experimental observation of cessation ofexpression dynamics before the mesoderm is segmentedinto somites [73]. This eventually gives way to a het-erogeneous steady state characterized by adjacent cellshaving alternating high and low clock gene expressions.The gap between these high and low values increases withtime to eventually produce a distinctive pattern that re-semble the stripes that arise due to polarization of eachsomite into rostral and caudal halves [as seen for large t in both Figs. 4 (a) and (b)]. In this asymptotic steadystate, the separation between the high and low valuesfor clock gene expression is greater than the amplitudeof the oscillations seen at lower values of t . Thus, wecan reproduce the entire sequence of dynamical transi-tions observed in the PSM during somitogenesis througha model incorporating an array of oscillators that inter-act via Notch-Delta signaling while “moving up” a mor-phogen gradient. IV. DISCUSSION
Somitogenesis is seen across all vertebrates, and recentevidence implies that mechanisms underlying it couldhave analogues even in segmentation of invertebrates,such as arthropods [5, 74]. It would appear that thereis an invariant set of mechanisms responsible for thisprocess, that differ only in terms of the specific identi-ties of the contributing molecular players across species.Thus, somitogenesis would in general involve (i) a cellular“clock”, (ii) means by which neighboring clocks commu-nicate, and (iii) a spatial gradient of signaling molecules,which introduces heterogeneity in the interactions be-tween the clocks. We have shown here that incorporatingthese three elements in a model of a PSM, that growsthrough the addition of cells at the posterior, reproducesthe sequence of invariant dynamical transitions seen insomitogenesis. These include the phenomenon of somitepolarization, provided that the genes determining thefates of cells constituting the rostral and caudal halvesof the mature somites [e.g., mesp-b in zebrafish [75–77]]are located downstream of the clock genes.While the roles of interacting clocks and that of mor- phogen gradients have been investigated individually inearlier studies, we provide here a framework to under-stand how these two work in tandem to give rise to thekey features associated with somitogenesis. In particular,our results shed light on the significance of the steepnessof the morphogen gradient. For instance, we may con-sider the consequences of a reduction in the steepnessleading to a linear profile for the morphogen gradientwhich can arise, for example, when the degradation rateis negligible. On replacing the exponential morphogengradient in our model with a linear one, we observe a verylong-lived transient state before the system converges toan inhomogeneous steady state. This therefore suggeststhat exponential gradients allow relatively rapid switch-ing between qualitatively distinct dynamical regimes [seeSI for results with linear gradient]. Hence, by varying thesteepness of the RA gradient experimentally it should bepossible to determine how the time required for matu-ration changes as a consequence. This is especially truein the case of the time interval between cessation of os-cillations and the polarization of the somites. As inter-cellular coupling is also known to regulate the period ofthe segmentation clock [78], it is possible that introduc-ing other morphogen gradients (such as, Wnt and FGF),that influence the strength of the coupling, can explainvariations in the rate at which somites form over time.Furthermore, the core assumption of our model, namelythat Notch-Delta coupling plays a crucial role in regulat-ing somitogenesis in the presence of a morphogen gradi-ent can be probed in experimental systems where Notchsignaling has been arrested.
ACKNOWLEDGMENTS
We would like to thank Krishnan Iyer, Jose Negrete Jr,Shubha Tole and Vikas Trivedi for valuable suggestions.SNM has been supported by the IMSc Complex Sys-tems Project (12 th Plan), and the Center of Excellencein Complex Systems and Data Science, both funded bythe Department of Atomic Energy, Government of India.The simulations and computations required for this workwere supported by High Performance Computing facility(Nandadevi and Satpura) of The Institute of Mathemat-ical Sciences, which is partially funded by DST. [1] M.-L. Dequ´eant and O. Pourqui´e, Nat Rev Genet , 370(2008).[2] C. Gomez, E. M. ¨Ozbudak, J. Wunderlich, D. Baumann,J. Lewis, and O. Pourqui´e, Nature , 335 (2008).[3] F. J. Vonk and M. K. Richardson, Nature , 282(2008).[4] S. Gilbert, Developmental Biology (Sinauer, 2013).[5] A. Stollewerk, M. Schoppmeier, and W. G. Damen, Na-ture , 863 (2003).[6] A. F. Sarrazin, A. D. Peel, and M. Averof, Science , 338 (2012).[7] O. Pourqui´e and P. P. Tam, Dev Cell , 619 (2001).[8] J. Cooke and E. C. Zeeman, J Theor Biol , 455 (1976).[9] R. E. Baker, S. Schnell, and P. Maini, Dev Biol , 116(2006).[10] M. Santill´an and M. C. Mackey, PLOS ONE , e1561(2008).[11] H. Nagahara, Y. Ma, Y. Takenaka, R. Kageyama, andK. Yoshikawa, Phys Rev E , 021906 (2009).[12] S. D. Hester, J. M. Belmonte, J. S. Gens, S. G. Clende- non, and J. A. Glazier, PLOS Comput Biol , e1002155(2011).[13] S. Ares, L. G. Morelli, D. J. J¨org, A. C. Oates, andF. J¨ulicher, Phys Rev Lett , 204101 (2012).[14] P. J. Murray, P. K. Maini, and R. E. Baker, Dev Biol , 407 (2013).[15] D. J. Jrg, L. G. Morelli, S. Ares, and F. Jlicher, PhysRev Lett , 174101 (2014).[16] G. Wiedermann, R. A. Bone, J. C. Silva, M. Bjorklund,P. J. Murray, and J. K. Dale, eLife , e05842 (2015).[17] O. Pourqui´e, Science , 328 (2003).[18] A. C. Oates, L. G. Morelli, and S. Ares, Development , 625 (2012).[19] I. Palmeirim, D. Henrique, D. Ish-Horowicz, andO. Pourqui´e, Cell , 639 (1997).[20] K. J. Dale and O. Pourqui´e, Bioessays , 72 (2000).[21] Y. Saga and H. Takeda, Nat Rev Genet , 835 (2001).[22] M. Maroto, J. K. Dale, M.-L. Dequeant, A.-C. Petit, andO. Pourquie, Int J Dev Biol , 309 (2003).[23] Y. Masamizu, T. Ohtsuka, Y. Takashima, H. Naga-hara, Y. Takenaka, K. Yoshikawa, H. Okamura, andR. Kageyama, Proc Natl Acad Sci USA , 1313 (2006).[24] I. H. Riedel-Kruse, C. M¨uller, and A. C. Oates, Science , 1911 (2007).[25] C. Schr¨oter, S. Ares, L. G. Morelli, A. Isakova, K. Hens,D. Soroldoni, M. Gajewski, F. J¨ulicher, S. J. Maerkl,B. Deplancke, et al. , PLOS Biol , e1001364 (2012).[26] A. B. Webb, I. M. Lengyel, D. J. J¨org, G. Valentin,F. J¨ulicher, L. G. Morelli, and A. C. Oates, eLife ,e08438 (2016).[27] Y.-J. Jiang, L. Smithers, and J. Lewis, Curr Biol , R868(1998).[28] Z. Ferjentsik, S. Hayashi, J. K. Dale, Y. Bessho, A. Her-reman, B. De Strooper, G. del Monte, J. L. de la Pompa,and M. Maroto, PLoS Genet , e1000662 (2009).[29] A. Hubaud and O. Pourqui, Nat Rev Mol Cell Bio ,709 (2014).[30] R. A. Conlon, A. G. Reaume, and J. Rossant, Develop-ment , 1533 (1995).[31] O. Pourqui´e, Curr Opin Genet Dev , 559 (1999).[32] J. Yun-Jin, B. L. Aerne, L. Smithers, C. Haddon, et al. ,Nature , 475 (2000).[33] E. C. Lai, Development , 965 (2004).[34] S. S. Huppert, M. X. G. Ilagan, B. De Strooper, andR. Kopan, Dev Cell , 677 (2005).[35] A. Mara and S. A. Holley, Trends Cell Biol , 593(2007).[36] R. Kageyama, Y. Masamizu, and Y. Niwa, Dev Dynam , 1403 (2007).[37] D. Sprinzak, A. Lakhanpal, L. LeBon, L. A. Santat, M. E.Fontes, G. A. Anderson, J. Garcia-Ojalvo, and M. B.Elowitz, Nature , 86 (2010).[38] D. Sprinzak, A. Lakhanpal, L. LeBon, J. Garcia-Ojalvo,and M. B. Elowitz, PLOS Comput Biol , e1002069(2011).[39] J. Lewis, Curr Biol , 1398 (2003).[40] F. Giudicelli and J. Lewis, Curr Opin Genet Dev , 0(2004).[41] K. Horikawa, K. Ishimatsu, E. Yoshimoto, S. Kondo, andH. Takeda, Nature , 719 (2006).[42] H. B. Tiedemann, E. Schneltzer, S. Zeiser, B. Hoesel,J. Beckers, G. K. Przemeck, and M. H. De Angelis, PLOSComput Biol , e1002586 (2012).[43] H. B. Tiedemann, E. Schneltzer, S. Zeiser, W. Wurst, J. Beckers, G. K. H. Przemeck, M. Hrab de Angelis, andD. Thieffry, PLOS Comput Biol , e1003843 (2014).[44] M. J. McGrew and O. Pourqui´e, Curr Opin Genet Dev , 487 (1998).[45] S. Gibb, M. Maroto, and J. K. Dale, Trends Cell Biol , 593 (2010).[46] J. Dubrulle, M. J. McGrew, and O. Pourqui´e, Cell ,219 (2001).[47] J. Dubrulle and O. Pourqui, Nature , 419 (2004).[48] J. Dubrulle and O. Pourqui´e, Development , 5783(2004).[49] J. Vermot and O. Pourqui´e, Nature , 215 (2005).[50] A. Aulehla, Gene Dev , 2060 (2004).[51] A. Aulehla, W. Wiegraebe, V. Baubet, M. B. Wahl,C. Deng, M. Taketo, M. Lewandoski, and O. Pourqui´e,Nat Cell Biol , 186 (2008).[52] L. Bajard, L. G. Morelli, S. Ares, J. Pecreaux, F. Julicher,and A. C. Oates, Development , 1381 (2014).[53] A. Goldbeter, D. Gonze, and O. Pourqui´e, Dev Dynam , 1495 (2007).[54] A. Goldbeter and O. Pourqui´e, J Theor Biol , 574(2008).[55] K. Mazzitello, C. Arizmendi, and H. Hentschel, PhysRev E , 021906 (2008).[56] D. J. J¨org, A. C. Oates, and F. J¨ulicher, Phys Biol ,05LT03 (2016).[57] A. D. Lander, Cell , 955 (2011).[58] R. Guantes and J. F. Poyatos, PLOS Comput Biol , e30(2006).[59] C. Takke and J. A. Campos-Ortega, Development ,3005 (1999).[60] A. C. Oates and R. K. Ho, Development , 2929(2002).[61] M. Maroto, J. K. Dale, M.-L. Dequeant, A.-C. Petit, andO. Pourquie, Int J Dev Biol , 309 (2005).[62] J. Gurdon and P.-Y. Bourillot, Nature , 797 (2001).[63] A. Aulehla and O. Pourqui´e, Cold Spring Harbor per-spectives in biology , a000869 (2010).[64] A. Sawada, M. Shinya, Y.-J. Jiang, A. Kawakami,A. Kuroiwa, and H. Takeda, Development , 4873(2001).[65] T. A. Moreno and C. Kintner, Dev Cell , 205 (2004).[66] A. D. Lander, Q. Nie, and F. Y. Wan, Dev Cell , 785(2002).[67] S. Bergmann, O. Sandler, H. Sberro, S. Shnider,E. Schejter, B.-Z. Shilo, and N. Barkai, PLOS Biol , 1(2007).[68] N. Barkai and B.-Z. Shilo, Cold Spring Harbor Perspec-tives in Biology (2009), 10.1101/cshperspect.a001990.[69] N. P. Shih, P. Fran¸cois, E. A. Delaune, and S. L.Amacher, Development , 1785 (2015).[70] R. D. del Corral, I. Olivera-Martinez, A. Goriely, E. Gale,M. Maden, and K. Storey, Neuron , 65 (2003).[71] R. D. del Corral and K. G. Storey, Bioessays , 857(2004).[72] M. Rhinn and P. Doll´e, Development , 843 (2012).[73] G. Tenin, D. Wright, Z. Ferjentsik, R. Bone, M. J. Mc-Grew, and M. Maroto, BMC Dev Biol , 24 (2010).[74] E. Clark, A. D. Peel, and M. Akam, Development (2019), 10.1242/dev.170480.[75] A. Sawada, A. Fritz, Y. Jiang, A. Yamamoto, K. Yamasu,A. Kuroiwa, Y. Saga, and H. Takeda, Development ,1691 (2000).[76] A. Nomura-Kitabayashi, Y. Takahashi, S. Kitajima, T. Inoue, H. Takeda, and Y. Saga, Development ,2473 (2002).[77] C.-Y. K. Kuan, D. Tannahill, G. M. Cook, and R. J.Keynes, Mech Develop , 1055 (2004).[78] L. Herrgen, S. Ares, L. G. Morelli, C. Schr¨oter,F. J¨ulicher, and A. C. Oates, Curr Biol , 1244 (2010). SUPPLEMENTARY INFORMATION
LIST OF SUPPLEMENTARY FIGURES
1. Fig S1: Schematic diagram of the genetic oscillator model and phase portrait of its dyanmics.2. Fig S2: Representative patterns of collective dynamics for a pair of coupled oscillators.3. Fig S3: Flow chart of the algorithm used to classify the collective dynamical patterns.4. Fig S4: Collective dynamics of a cellular array represented in terms of the bound Notch concentration ( N b )5. Fig S5: Temporal evolution of the concentration of a morphogen along the gradient as sensed by a particularcell in the growing PSM.6. Fig S6: Temporal evolution of the dynamical variables in a pair of coupled oscillators responding to an expo-nential gradient of morphogen concentration.7. Fig S7: Effect of delay in Notch-Delta mediated interaction on the collective dynamics of a pair of coupledoscillators responding to an exponential gradient of morphogen concentration cell signaling. MODELING THE DYNAMICS OF CLOCK GENE EXPRESSION IN CELLS COUPLED VIANOTCH-DELTA SIGNALING
In the model presented in the main text, we consider a two-component genetic oscillator [58], comprising anactivator gene x and an inhibitor gene y . The parameter values of the model (see Table S1) have been chosen suchthat it exhibits autonomous oscillatory activity. The interactions between these two genes are schematically shownin Fig. S1 [left]. Each uncoupled oscillator consists of two variables X and Y that represent the concentrations of theproducts expressed by genes x and y , respectively. The trajectory of an uncoupled oscillator in the X - Y phase planeis shown in Fig. S1 [right], along with the nullclines ˙ X = 0 and ˙ Y = 0. Fig. S1 [left] also displays the nature of theinteractions between genes x and y , where the variables N b , N and D describe the Notch-Delta coupling betweencells [37].In a system of coupled oscillators, we observe a variety of synchronization behavior (examples of which are shownin Fig. S2). These can be categorized into 6 principal collective dynamical patterns, namely Inhomogeneous SteadyState (ISS), Exact Synchronization (ES), Homogeneous Steady State (HSS), Inhomogeneous In-phase Synchronization(IIS), Anti-Phase Synchronization (APS) and Other synchronization patterns (OTH). In order to classify all observeddynamical patterns into these 6 classes, we use a set of order parameters, as illustrated in the flow chart shown inFig. S3. We have verified that the results are robust with respect to small changes in the thresholds used to determinethese patterns.We would like to point out that, while in the main text we have exclusively used the variable Y to illustrate thedynamical transitions, qualitatively identical behavior can be seen using other variables, such as N b (as can be seenon comparing Fig. S4 with Fig. 4 of the main text). parameter a b c e β N γ k µ value 16.0 200.0 20.0 10.0 5.0 1.0 1.0 1.0TABLE S1. Model parameter values used for all simulations in the main text (unless specified otherwise). FIG. S1. [left] Schematic diagram of the interactions resulting in gene expression oscillations in the model used in the main text.The system comprises an activator gene x and an inhibitor gene y that yield the gene products X and Y , respectively. Thesein turn inhibit the expression of Delta gene d , resulting in suppression of the production of Delta ligands D . Each cell containsNotch receptors N that, when bound with Delta ligands from another cell, causes the Notch Intracellular Domain (NICD,represented through the proxy variable N b ) to cleave off and act as transcriptional factors that upregulate the expressionof the downstream genes x and y . [right] The dynamics of expression levels X and Y for the activator and inhibitor genes,respectively, in a single uncoupled oscillator, shown in terms of the trajectory of the limit cycle (broken curve) and the nullclines(red: ˙ X = 0, blue: ˙ Y = 0) obtained from Eqs. (1)-(2) in the main text, with the parameter values shown in Table S1. t Y ISS t Y HSS t Y APS t Y ES t Y IIS t Y OTH
FIG. S2. Typical time series of the inhibitor gene expression Y for a pair of cells coupled to each other (each exhibitingoscillations in isolation, as shown in Fig. S1). Different dynamical regimes obtained by varying the coupling parameters f , g and h that are shown here correspond to (top row, left) Inhomogeneous Steady State (ISS), (bottom row, left) ExactSynchronization (ES), (top row, center) Homogeneous Steady State (HSS), (bottom row, center) Inhomogeneous In-phaseSynchronization (IIS), (top row, right) Anti-Phase Synchronization (APS) and (bottom row, right) Other synchronizationpatterns (OTH). The activity of the two cells are represented using two different colors in each panel. Note that the two curvesoverlap in the ES and HSS regimes. (cid:104) σ t ( Y i ) (cid:105) i (cid:104) Y i Y j (cid:105) t ES > . APS < − . IIS . ↔ . Others − . ↔ . > − σ i ( (cid:104) Y i (cid:105) t ) HSS < − ISS > − < − FIG. S3. Flow chart illustrating the algorithm used to classify the collective dynamical patterns obtained for a system oftwo coupled oscillators. To distinguish between oscillating and non-oscillating patterns, we use (cid:104) σ t ( Y i ) (cid:105) i which is the meantemporal variance of the time series of Y , calculated over the two oscillators. To determine whether the steady states that theoscillators have reached are the same (corresponding to HSS) or different (corresponding to ISS), we compute the variance ofthe mean values for the two time series, σ i ( (cid:104) Y i (cid:105) t ). To distinguish between the oscillating patterns, we use the equal time linearcorrelation between the two time series, (cid:104) Y i Y j (cid:105) t . The classification is robust with respect to small changes in the values of thethresholds, which are displayed in the figure. c e ll i nde x a b b t (a.u.) b FIG. S4. Collective dynamics of a cellular array responding to an exponential gradient of morphogen concentration reproducesthe spatio-temporal evolution of PSM activity seen during somitogenesis, similar to Fig. 4 in the main text, but showing boundNotch ( N b ) instead of the inhibitor variable Y . MODELING THE MORPHOGEN GRADIENT IN A GROWING PSM
A crucial element of our model is the gradient of morphogen concentration across the presomitic mesoderm (PSM).We have explicitly considered the morphogen to be Retinoic acid (RA) whose concentration is highest at the anteriorend of the PSM. Experimental evidence suggests that as the anteriorly located newly formed somites mature, each ofthem serve as a source of RA [70–72]. This suggests that a source, diffusion and linear degradation (SDD) model forRA will provide an appropriate description for the variation of the morphogen concentration across space and time.In particular, the source grows in strength as the number of mature somites increases over time. If one considers thePSM to be growing in discrete steps, with each new mature somite (which secretes RA at a rate φ ) being added aftera time interval τ d , the concentration C i of the morphogen at each putative somite in position i on the anteroposterioraxis of the growing PSM can be described by the differential equation: dC i dt = R ( t ) δ i, + D ( C i +1 + C i − − C i ) − C i τ m , where the time-varying production term R ( t ) increases in a step-like manner by φ after each interval of duration τ d , δ i, is a Kronecker delta function indicating that the source is at the anterior end of the domain (using the simplifyingassumption of a point source rather than a distributed one), and D and τ m are the effective diffusion rate and averagelifetime of the RA molecules, respectively. We note that qualitatively identical results are obtained if, instead ofincreasing in discrete time steps, the production rate R ( t ) changes in a continuous fashion, for example, as describedby the differential equation dRdt = φτ d . Fig. S5 shows the resulting exponential profile of the morphogen concentration that will be experienced by a particularcell in the expanding PSM, which validates the use of an exponential profile for the morphogen in the main text.
FIG. S5. Concentration of the morphogen gradient as sensed by a particular cell in the growing PSM with time, shown for (a)different values of φ for τ m = 10, and (b) different values of τ m for φ = 0 .
4. The simulation domain comprises N = 10 putativesomites at any given time, while the other parameter values are D = 1, τ d = 5 and R ( t = 0) = 1 . THE EFFECT OF DELAY ON THE DYANMICS OF COUPLED OSCILLATORS IN THE PRESENCE OFA MORPHOGEN GRADIENT
In the main text, we have described the dynamics of clock gene expressions in cells interacting via Notch-Deltasignaling (modulated by a morphogen gradient) where communication delay between the cells is assumed to beinsignificant. In Fig. 3 of the main text, we show the dynamical transitions in a pair of coupled cells by focusing onthe inhibitor gene product concentration Y . We would like to point out that similar dynamics is observed in the othersystem variables (see Fig. S6). The only exception is the free Notch receptor concentration N which, unlike the othervariables, grows in an unregulated manner.In order to investigate the effect of an explicit delay in the system, we incorporate a delay of duration τ d in both(i) the regulation of X and Y by N b and (ii) the repression of the ligand Delta by X and Y . In Fig. S7, we showthe effect of increasing delay, which is expressed in terms of the relative magnitude τ d /τ P , where τ P is the period ofthe uncoupled oscillator, on the dynamics of the system. We observe the same qualitative nature for the dynamicaltransitions as seen in Fig. 3 (d) of the main text, suggesting that our results are robust with respect to incorporationof signaling delays. (a) (b) (c) b (d)
100 200 30001 t (a.u.)D (e)