Hopf bifurcation in the presomitic mesoderm during the mouse segmentation

ARTICLE IN PRESS Journal of Theoretical Biology 259 (2009) 176–189

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Hopf bifurcation in the presomitic mesoderm during the mouse segmentation Aitor Gonza´lez a,b,, Ryoichiro Kageyama a,b a b

Institute for Virus Research, Kyoto University, Japan Japan Science and Technology Agency, CREST, Shogoin-Kawahara, Sakyo-ku Kyoto 606-8507, Japan

a r t i c l e in f o

a b s t r a c t

Article history: Received 10 February 2009 Accepted 11 February 2009 Available online 21 February 2009

Vertebrae and ribs arise from embryonic tissues called somites. Somites arise sequentially from the unsegmented embryo tail, called presomitic mesoderm (PSM). The pace of somite formation is controlled by gene products such as hairy and enhancer of split 7 (Hes7) whose expression oscillates in the PSM. In addition to the cyclic genes, there is a gradient of fibroblast growth factor 8 (Fgf8) mRNA from posterior to anterior PSM. Recent experiments have shown that in the absence of Fgf signaling, Hes7 oscillations in the anterior and posterior PSM are lost. On the other hand, Notch mutants reduce the amplitude of posterior Hes7 oscillations and abolish anterior Hes7 oscillations. To understand these phenotypes, we delineated and simulated a logical and a delay differential equation (DDE) model with similar network topology in wild-type and mutant situations. Both models reproduced most wild-type and mutant phenotypes suggesting that the chosen topology is robust to explain these phenotypes. Numerical continuation of the model showed that even in the wild-type situation, the system changed from sustained to damped, i.e. a Hopf bifurcation occurred, when the Fgf concentration decreased in the PSM. This numerical continuation analysis further indicated that the most sensitive parameters for the oscillations are the parameters of Hes7 followed by those of Lunatic fringe (Lfng) and Notch1. In the wild-type, the damping of Hes7 oscillations was not so strong so that cells reached the new somites before they lose Hes7 oscillations. By contrast, in the fibroblast growth factor receptor 1 (Fgfr1) conditional knock-out (cKO) mutant simulation, Notch signaling was not able to maintain sustained Hes7 oscillations. Our analysis suggests that Fgf signaling makes cells enter an oscillatory state of Hes7 expression. After moving to the anterior PSM, where Fgf signaling is missing, Notch signaling compensates the damping of Hes7 oscillations in the anterior PSM. & 2009 Elsevier Ltd. All rights reserved.

Keywords: Dynamical modeling Mouse segmentation Notch signaling Fgf signaling bHLH gene oscillations

1. Introduction Vertebrae and ribs arise from embryonic tissues called somites. In mouse embryos, somites appear as blocks from the unsegmented tail, which is called presomitic mesoderm (PSM). Cells do not significantly move in the PSM relative to each other. However, because cells divide mainly in the tail bud and differentiate to somites in the anterior PSM, a PSM cell can be thought as moving relative to the tail bud towards the newest somite (Fig. 1(a)) (reviewed by Gonza´lez and Kageyama, 2007). Periodic somite formation requires cyclic expression in the PSM of the so-called cyclic genes. Microarray and in situ hybridization data have shown that many cyclic genes exist (Deque´ant et al., 2006). Many of those cyclic genes belong to the Notch, fibroblast growth factor (Fgf) and Wnt signaling pathways.

 Corresponding author at: Institute for Virus Research, Kyoto University, Japan. Tel.: +8175 751 4013. E-mail address: [email protected] (A. Gonza´lez).

0022-5193/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2009.02.007

It is currently not clear whether a single master clock or several interlocked clocks drive the oscillations. An important candidate for master clock is the Notch pathway gene hairy and enhancer of split 7 (Hes7). In addition of being required for the correct mammalian segmentation, the Hes7 protein exerts a negative feedback on its own expression (Fig. 1(b)) (Bessho et al., 2003). Within the Notch signaling pathway, there are other important feedbacks like that where Hes7 inhibits transcriptionally the glycosyltransferase Lunatic fringe (Lfng). Lfng in turn glycosylates the Notch receptor, which upon activation induces Hes7 and Lfng expression (Bessho et al., 2003; Morimoto et al., 2005). In addition to the cyclic genes, there are three cellular pathways activated as a gradient along the PSM, the Fgf and Wnt pathways (peaks in the tail bud), and the retinoic acid (RA) pathway (peak in the new somite) (reviewed by Gridley, 2006). In the most broadly accepted working model, the clock-andwavefront model, those gradients are thought to stop the expression oscillations of the cyclic genes and so, to define a new somite as a function of the oscillation phase of the cyclic genes (Cooke and Zeeman, 1976). In recent years, several links

ARTICLE IN PRESS ´lez, R. Kageyama / Journal of Theoretical Biology 259 (2009) 176–189 A. Gonza

177

Fig. 1. (a) PSM of an E10.5 mouse embryo stained with Hes7 intronic probe (blue), and (b) regulatory components modeled in this paper (evidence and references of the Notch module are in Table 1). (a) In the PSM, there is a gradient of Fgf8 expression decreasing from the tail bud to the forming somite. The PSM grows posteriorly and segmentates anteriorly, so that cells change their position in the anterior direction relative to the tail bud and the active Fgfr1 gradient peak. (b) The Fgf/RA module was modeled as Goldbeter et al. (2007). Retinoic acid (RA) is catalyzed by the enzyme Raldh2 and hydrolyzed by the enzyme Cyp26. Fgf8 is a secreted molecule that activates the Fgfr1 receptor, which then induces the expression of Cyp26 gene. On the other hand, RA represses the expression of the Fgf8 mRNA (Goldbeter et al., 2007). In the Notch module, there are three components. The bHLH gene Hes7 is induced by the Fgfr1 and Notch1 receptors (references in Table 1). Hes7 binds and inhibits its own promoter as well as the Lfng promoter. Lfng is a glycosyltransferase enzyme that glycosylates the Notch1 receptor. This results in the inhibition of the Notch1 receptor, although there are contrasting views about this aspect (see Discussion). The Notch1 receptor is activated by a Dll1 ligand, which is present throughout the PSM, so that we have omitted it from our model. The Notch module is formalized in a similar way as Rodriguez-Gonzalez et al. (2007) (Eqs. (2)). Normal and blunt arrows stand for activating and inhibiting interactions, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

between the gradients and cyclic genes were found. For instance, Aulehla et al. (2003) reported that Wnt signaling is necessary for the cyclic gene Axin2 expression that in turn represses Wnt signaling. Wahl et al. (2007) reported that the fibroblast growth factor receptor 1 (Fgfr1) is required for expression of several cyclic genes of the Notch and Wnt signaling pathway. Recently, a new link between the Hes7 molecular clock and the Fgf8 wavefront was found. It was found that Fgf is required for Hes7 oscillations (Niwa et al., 2007). Interestingly, this requirement was not uniform along the PSM and also not in the duration. Short treatment (2 h) of the PSM with Fgf inhibitors abolished posterior but not anterior Hes7 oscillations. By contrast, long treatment (6 h) of the PSM with Fgf inhibitors and Fgfr1 conditional knock-out (cKO) embryos extended this perturbation of Hes7 oscillations to the anterior PSM. Finally, short and long treatments of the PSM with the Notch inhibitor DAPT, and genetic perturbation in Dll1 knock-out (KO) (ligand of Notch1 receptor), Rbpj cKO (transcriptional effector of the Notch pathway) and Lfng transgenic embryos affected Hes7 oscillations only in the anterior PSM (Serth et al., 2003; Niwa et al., 2007). To try to understand the links between the PSM gradients, different theoretical models have been presented. In one of those models, the Wnt pathway regulated cyclic genes belonging to the Wnt pathway as well as cyclic genes of the Notch pathway (Rodriguez-Gonzalez et al., 2007; Santilla´n and Mackey, 2008). In another model by Goldbeter and Pourquie´ (2008), the Fgf pathway activated Axin2 and Dusp6, while the Wnt pathway regulated both Axin2 and the Notch intracellular domain. In another model, the Fgf signaling pathway increased the degradation rate of the Hes7 mRNA, so that the oscillation period became longer as the cells moved towards the anterior PSM (Tiedemann et al., 2007). In a recent model of the segmentation in zebrafish, parameters that increase the period of the oscillations and thus might be regulated by the PSM gradient were shown (Uriu et al., 2009). However, those models were not challenged by simulations of loss-of-function mutations. Furthermore, all of these models were based on differential equations that required a great detail about parameter values, which is not available.

To overcome these limitations, we proposed a model of the crucial molecules in the intracellular Notch pathway including Hes7, Lfng and Notch1 using both logical and delay differential equations (DDEs). We used this model to formalize the observations of Hes7 phenotypes in Fgf and Notch mutants. The model was evaluated by simulating wild-type and mutant situations.

2. Materials and methods 2.1. Simplifications and limitations of the model Microarray data have found many cyclic genes that belong to at least one of the three signaling pathways Wnt, Fgf and Notch (Deque´ant et al., 2006). Also many negative feedback loops might drive the oscillation of cyclic genes (Gonza´lez and Kageyama, 2007). It is not the aim of the present model to integrate all interactions as in the case of previous works (Goldbeter and Pourquie´, 2008). Here, we wanted to test whether the Notch subnetwork activated by the Fgf signaling is sufficient to explain the observations made by Niwa et al. (2007). Though the ligands of Notch1, i.e. Dll1 and Dll3, are required for the correct segmentation, their expression seems to be ubiquitous or at least non-oscillatory in the PSM (del Barco Barrantes et al., 1999; Niwa et al., 2007; Shifley et al., 2008). Therefore, we assumed that Dll1 is implicitly present all along the PSM, and the basal level of Notch1 was active (in the absence of Lfng). For that reason, and because the model only included a single cell moving along the PSM, we did not investigate the mechanism of cell synchronization. The Fgf signaling pathway was not modeled in detail to simplify the analysis of mutants. The regulation of Hes7 by Fgf was assumed to be at the transcriptional level, because this is the simplest explanation for the loss of Hes7 intronic mRNA expression in the absence of Fgf signaling (Niwa et al., 2007). We did not take into account separate variables for the mRNA and protein in the models presented here in order to reduce the number of variables and parameters. We only took into account a single variable that could be seen as the activity of the gene.

ARTICLE IN PRESS 178

´lez, R. Kageyama / Journal of Theoretical Biology 259 (2009) 176–189 A. Gonza

2.2. Delineation of the model interactions Fgfr1 was assumed to activate the synthesis of Hes7 as this is the simplest explanation for the observation that Fgfr1 loss-offunction mutation resulted in loss of Hes7 mRNA expression (Niwa et al., 2007). We simplified this interaction by drawing a positive interaction from Fgfr1 to Hes7. On the other hand, it is known that Notch signaling regulates Hes7 at the transcriptional level (Bessho et al., 2003). For the combined regulation of Hes7 by Fgf and Notch signaling, we specified a logical additive function because in the absence of Notch signaling, Hes7 mRNA is still expressed (Niwa et al., 2007). The inhibition of the Hes7 promoter by Hes7 protein was represented by a logical multiplicative function, because expression of Hes7 was able to block the upregulation of the Hes7 promoter by Notch signaling (Bessho et al., 2003). Notch1 was assumed to be constitutively upregulated in the absence of Lfng, as we assumed that Dll1 was present throughout the PSM (Section 2.1). Although Lfng is an enzyme, we represented the inhibition of Notch1 by Lfng through Boolean and Hill functions, because these functions well represented the fact that Notch1 is active in the absence and inactive in the presence of Lfng (Niwa et al., 2007). 2.3. Estimation of delays for Hes7, Lfng and Notch intracellular domain effects 2.3.1. Estimation of delays based on in situ hybridization data Fig. 5 presented by Bessho et al. (2003) showed expression patterns of intronic and mature Hes7 mRNA and Hes7 protein in three different phases of the oscillation pattern. In those data, Bessho et al. (2003) presented the proportion of embryos in each oscillation phase, i.e. 8/21 embryos in phase I, 5/21 in phase II and 8/21 in phase III. In average, ð8 þ 5Þ=2 ¼ 6:5=21 embryos were in the transition I/II, 6.5/21 in the transition II/III and 8/21 in the transition III/I. Assuming a segmentation period of 120 min, then the transitions I/II and II/III lasted 6:5=21  120 min ¼ 37 min and the transition III/I lasted 8=21  120 min ¼ 46 min. In phase I, the Hes7 mRNA intron but not the Hes7 protein were observed while in phase II, the Hes7 protein but not the Hes7 mRNA intron were observed. This indicated that the transcription/translation delay of Hes7 was around 37 min according to in situ hybridization data. Though in these data, the intronic pattern of Lfng was not shown, its mRNA pattern looked similar, so that we assumed the transcription/translation delay of Lfng to be similar to that of Hes7. 2.3.2. Estimation of delays based on the sequence The mouse Hes7 gene has 2808 base pairs (bp), 3 introns and 225 amino acids (aa), whereas the Lfng gene has 2282 bp, 7 introns and 378 aa. We assumed a transcription rate of 1200 bp/min (O’Brien and Lis, 1993), which results in a transcriptional delay of 2.34 min for Hes7 and 1.9 min for Lfng. The splicing was assumed to have a half-life of 0.4–7.5 min/intron (Audibert et al., 2002), which results in a splicing delay of 1.2–22.5 min for Hes7 (3 introns) and 2.8–52.5 min for Lfng (7 introns). Nuclear export might take around 4 min (Audibert et al., 2002). Translation takes around 3.5 aa/s (210 aa/min) (Ujvari et al., 2001), so that the translation delay is 1.1 min for Hes7 and 1.8 min for Lfng. Finally, Hes7 needs to return to the nucleus to exert its function, so that a further nuclear import delay of around 3–5 min is taken into account (Zhu et al., 2007). In summary, based on sequence data, the transcription/translation delay of Hes7 is in the range 11.64–37 min, and that of Lfng in the range 10.5–60.2 min. According to the calculations based on the in situ hybridization data, it seems that there might be a transcription/translation delay of around 37 min for both Hes7 and Lfng. In cell culture, the delay between mature mRNA and protein is around 20 min,

so that the actual transcription/translation delay should be larger (Hirata et al., 2002). For the sake of comparison, in zebrafish in situ hybridization and immunostaining data showed that the transcription/translation delay is about 30 min for DeltaC (Giudicelli et al., 2007). 2.4. Normalization We started from these initial equations:   krh h dh nrn ln2  h, ¼ ah m Fgfr1 þ ð1  mÞ rn dt kn þ nrn krh h þ hrT h th h

rh

kh dl nrn ln2  l, ¼ al rn dt kn þ nrn krh h þ hrT hh tl r

k l dn ln2 n, ¼ an rl l rl  tn dt kl þ lT l

where h, l, n represented Hes7, Lfng and Notch1 activities, respectively. The delayed variables hT h ; lT l stand for hðt  T h Þ; lðt  T l Þ, respectively. The parameters were the synthesis rates ah ; al ; an , the half-lives th ; tl ; tn , the relative contribution m of Fgf signaling to the Hes7 promoter and the threshold values kh ; kl ; kn (equivalent to the Michaelis–Menten coefficient). We defined the normalized variables: h~ ¼ ln 2=ðah th Þh; ~l ¼ ln 2=ðal tl Þl; n~ ¼ ln 2=ðan tn Þn. We also defined a normalized yh ¼ kh ln 2=ðah th Þ; yl ¼ kl ln 2=ðal tl Þ; yn ¼ kn ln 2= threshold ðan tn Þ. After substituting in the previous equations, rearranging a bit and throwing the tildes, we obtained Eqs. (2). 2.5. Model of the Fgf and RA pathways We used a model from the literature to represent the Fgf/RA module (Goldbeter et al., 2007). We coupled this model to Eqs. (2) with minor notational changes: dRA ¼ vs1  kd1 C RA  kd5 RA, dt dMC F2 ¼ V 0 þ V sC 2  K d3 M C , dt KA þ F2 dC ¼ ks2 MC  kd2 C, ct dF K2 ¼ ks3 M F 2 I 2  kd4 F, dt K I þ RA where RA, MC , C and F variables correspond to the RA, Cyc26 mRNA, Cyc26 protein and Fgf8 protein concentrations/levels, respectively. The RA synthesis rate is given by vs1 ¼ ks1 RALDH20 ð1  x=LPSM Þ. The position in the PSM is given by x ¼ LPSM ð1  t= PPSM Þ. We assumed that the simulation time is equal to the time needed by a cell to cross the PSM (see below). We defined the asymmetric Fgf8 mRNA M f ¼ M 0 x=LPSM , the Fgf receptor occupancy Fgfr1 ¼ F=ðF þ K r2 Þ, and the RA receptor occupancy RAR ¼ RA=ðRA þ Kr 1 Þ. The parameter values were chosen from Goldbeter et al. (2007): ar0 ¼ 1; r 0 ¼ 7:1; M 0 ¼ 5; RALDH20 ¼ 7:1; ks1 ¼ 1; ks2 ¼ 1; ks3 ¼ 1; kd1 ¼ 1; kd2 ¼ 0:28; kd3 ¼ 1; kd4 ¼ 1; kd5 ¼ 0; V 0 ¼ 0:365; V sC ¼ 7:1; K r1 ¼ 1; K A ¼ 0:2; K I ¼ 0:3. 2.6. Estimation of the time required by a cell to move across the PSM Though cells do not significantly move relative to each other, the new somite approaches any PSM cell at the same rate as the

ARTICLE IN PRESS ´lez, R. Kageyama / Journal of Theoretical Biology 259 (2009) 176–189 A. Gonza

determination wavefront. At embryonic stage E10.5, the segmentation period is about 120 min and the new somites are about 145 mm in length (Tam, 1981). This implies that the determination wavefront moves at a rate of about 1.2 mm/min. At this stage, the length of the PSM is about 813 mm, so that a given cell will take around 683 min to reach the anterior border of the PSM (Tam, 1981). Based on these arguments, the simulation time of 683 min correlates with the PSM position and also with the Fgf signaling levels. For our parameter values, the period was around 130 min, so that we simulated the model for a little bit longer time, namely 720 min. 2.7. Software The logical model was defined and simulated with GINsim (Gonza´lez Gonza´lez et al., 2006). The DDE model was simulated with XPPAUT (Ermentrout, 2002). The numerical continuations were carried out with PDDECONT (Szalai, 2005). The results of the numerical continuations were confirmed with DDEBIFTOOL (Engelborghs et al., 2002) and by linear stability analysis (MacDonald, 1989) (data not shown).

the domain of the integer variables. The activities of Fgfr1, Hes7 and Lfng were associated with two values. On the other hand, the activity of Notch1 was associated with three values, a fully activated level 2, where Notch1 induces both Lfng and Hes7, and an intermediate level 1, where only Lfng but not Hes7 is activated, and an inactive level 0. This was justified by the observation that in the absence of Fgf signaling, Notch1 is still able to upregulate Lfng but not Hes7 (Niwa et al., 2007). The domain of the state variables was defined as F ¼ f0; 1g; H ¼ f0; 1g; L ¼ f0; 1g; N ¼ f0; 1; 2g, and the thresholds as: 0oyF o1; 0oyH o1; 0oyL o1; 0oyN;1 o1oyN;2 o2. The temporal evolution of the state variables at time t, i.e. F t ; Ht ; Lt ; N t was governed by logical parameters kF ; kH ; kL ; kN defined in the same domains as the corresponding state variables, i.e. kF ¼ f0; 1g; kH ¼ f0; 1g; kL ¼ f0; 1g; kN ¼ f0; 1; 2g, as a function of these Boolean equations: kF ¼ 0,

In situ hybridization in Fig. 1 was performed as previously described (Bessho et al., 2001).

3. Results

(

kH ðF t ; Ht ; Nt Þ ¼ ( kL ðHt ; N t Þ ¼ ( kN ðLt Þ ¼

2.8. In situ hybridization

179

1

if s ðHt ; yH Þ

0

else;

and

ðsþ ðF t ; yF Þ or sþ ðNt ; yN;2 ÞÞ ¼ true;

1

if s ðHt ; yH Þ and sþ ðN t ; yN;1 Þ ¼ true;

0

else;

2

if s ðLt ; yL Þ ¼ true;

1

else;

(1)

where the Boolean functions sþ ; s were defined as ( true X t 4yX ; sþ ðX t ; yX Þ ¼ false X t oyX ; s ðX t ; yX Þ ¼ not sþ ðX t ; yX Þ.

3.1. Logical multivalued model In developmental biology, evidence is mainly based on genetic experiments, i.e. genetic perturbation followed by assessment of either overexpression or lack of expression of another gene. This kind of causal relations can be most naturally formalized with logical multivalued models, where the levels of a gene are represented by multivalued logical variables with integer values 0; 1; 2; . . . ; and the effects of combinations of interactions by logical functions of the regulator values (Thomas and D’Ari, 1990; Thieffry and Thomas, 1995). To simplify the analysis of the logical model, we focused on the Notch module (yellow subgraph in Fig. 1(b)) under the control of Fgfr1. Therefore, we started our investigation by collecting genetic data and creating a logical multivalued model about the regulation of Hes7 (Table 1 and Section 2.2). For each gene/protein of the Notch module and Fgfr1, we defined integer variables F; H; L; N representing the state of Fgfr1, Hes7, Lfng and Notch1 relative to theoretical real-valued thresholds yF ; yH ; yL ; yN;1 ; yN;2 . Note that these thresholds were implicit in

Based on the previously defined logical parameters, the state of each variable X tþT X with X 2 fF; H; L; Ng at the future time t þ T X was computed as 8 X t þ 1 if kX  X t 40; > > < if kX  X t ¼ 0; X tþT X ¼ X t > > : X  1 if k  X o0: t

X

t

Note that kN could not take the value 0, because we assumed that Notch1 is able to reach an equilibrium with Lfng at value 1. At this equilibrium, the level 1 of Notch1 maintained Lfng expression, although Lfng did not allow Notch1 to reach its highest level 2. This was justified by the observation that ubiquitous Lfng expression was observed in the PSM of some mutants like Fgfr1 cKO. The delay implicit in transcription, translation, etc. was taken into account by the delays T F ; T H ; T L ; T N . These delays specify the commutation orders. The commutation order defines the order by which the variables change, when two or more variables are called to change. In the fully asynchronous simulation, all possible commutation orders are simulated. Note that this method does

Table 1 Regulatory components of the Notch module (Fig. 1) (1st column), regulatory interactions (2nd column) and experimental evidence (3rd column). Variable

Interaction

Fgfr1 Hes7

Lfng

Notch1

Hes7 a Hes7 Fgfr1!Hes7 Notch1!Hes7 Notch1!Lfng Hes7 a Lfng Lfng a Notch1

Experimental evidence The mRNA of the Fgf8 gene, the ligand of the Fgf receptor (Fgfr1), is distributed as a gradient from the posterior to the anterior PSM (Dubrulle and Pourquie´, 2004) Hes7-null embryos display homogeneous Hes7 transcription as assessed by a marker gene under the control of the Hes7 promoter. By contrast, Hes7 protein stabilized with a protease inhibitor abolishes Hes7 transcription (Bessho et al., 2003) Inhibition of Fgf signaling abolishes posterior Hes7 expression (Niwa et al., 2007) Inhibition of Notch signaling abolishes anterior Hes7 expression (Niwa et al., 2007) Inhibition of Notch signaling abolishes Lfng expression (Dale et al., 2006) Hes7-null embryos display homogeneous Lfng expression, whereas Hes7 protein stabilization abolishes Lfng expression (Bessho et al., 2003) In Lfng KO mouse, Notch1 activation is found throughout the PSM (Morimoto et al., 2005), whereas Lfng overexpression in chick results in loss of Notch targets (Dale et al., 2003)

ARTICLE IN PRESS 180

´lez, R. Kageyama / Journal of Theoretical Biology 259 (2009) 176–189 A. Gonza

not need numerical values for the delays. Because all possibilities were taken into account, two or more outgoing arcs (or successor states) were common for each state (Fig. 2). This model was formalized and simulated using the software GINsim (Gonza´lez Gonza´lez et al., 2006). To validate this logical model, we carried out wild-type and mutant simulations (Fig. 2). The simulation of this logical model takes the form of a socalled state transition graph, where nodes represent gene expression states and arcs represent transitions between gene expression states. More specifically, each node was composed of

four values that represented the expression/activation level of the variables F t ; Ht ; Lt ; N t at time t. In this graph, oscillations occurred, whenever there is a path from a state to itself. For instance, in the transitions 1011 ! 1111 ! 1011, the variable of Hes7 (second value) oscillated between values 0 and 1. To assess the validity of this model, we carried out simulations of the wild-type and mutant situations with a given initial state (square states in Fig. 2). To model the Fgfr1 gradient in the logical model, an initial state of F t ¼ 1 and a basal parameter value of kF ¼ 0 were set. By this way, the F t variable changed from value 1 to value 0 in a

Fig. 2. State transition graphs of wild-type (a) and mutant (b–e) situations for an initial state (square box) of the logical parameters in Eqs. (1). Each node was composed of four values that corresponded to the expression/activation levels of Fgfr1, Hes7, Lfng and Notch1. Nodes and arcs corresponded to system states and transitions between states, respectively. States were updated asynchronously, so that two or more outgoing arcs were observed. (a) In the wild-type simulation, Hes7 oscillated in the region with high (black font) and low (white font) Fgf levels. (b) Fgfr1 cKO simulation resulted in a complete inactive system. (c) Rbpj cKO simulation resulted in the inactivation of Hes7 oscillations only when the Fgf signaling levels were low (white font states). (d and e) Lfng KO simulation allowed the system to oscillate independently of Fgfr1 activation levels, whereas Hes7 KO simulation completely abolished system oscillations. KO mutant simulations were carried out by fixing the corresponding variables to value 0. Color codes: Green ¼ strong connected component (SCC), red ¼ stable states and gray ¼ neither SCC nor stable states. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

ARTICLE IN PRESS ´lez, R. Kageyama / Journal of Theoretical Biology 259 (2009) 176–189 A. Gonza

similar way as a given cells would sense decreasing levels of Fgf signaling when moving from the posterior to the anterior PSM. The wild-type simulation of the logical model showed that Hes7 oscillated for high ðFt ¼ 1Þ and low ðFt ¼ 0Þ Fgf signaling levels (Fig. 2(a)). Then, we simulated the perturbation experiments that we wanted to explain here. To simulate the Fgfr1 cKO mutant, we blocked the variable of Fgfr1 (first digit) to value 0. This simulation resulted in a single state 0011, where Hes7 is blocked at value 0 and Lfng and Notch1 at value 1 (Fig. 2(b)). Next, we simulated a Notch mutant by fixing the variable of Notch1 (last digit) to value 0. Also we carried out mutant simulations of the Lfng KO and Hes7 KO mutants by fixing the corresponding variables to value 0. Whereas the Lfng KO simulations showed oscillations (though somehow different from the wild-type), the Hes7 KO situation completely lacked oscillations in agreement with experimental observations (Figs. 2(d,e)). In summary, the logical model reproduced wild-type and mutant observations, including the PSM position specific phenotypes of the Fgf and Notch perturbation experiments. In the previous simulations, an initial state was chosen. To evaluate the robustness of our model to initial conditions, we simulated all possible states of the system, i.e. 23  3 ¼ 24 states for the wild-type (Fig. 3). To ease the graph visualization, we compressed the strong connected components (SCCs). The SCCs are subgraphs, where any state can be reached from any state, so that they include cyclic states. Again, we simulated wild-type and mutant simulations. These simulations also showed that the wildtype model oscillated in the presence of high and low Fgf levels (green SCCs in Fig. 3(a)). In the case of the Fgfr1 cKO mutant simulations, the system underwent some transient oscillations before reaching the steady state 0011 (SCC cc-0 in Fig. 3(b)). These oscillations were transient and would be eliminated if a presimulation was carried out. In the case of the Notch KO mutant, the system oscillated only in the presence but not in the absence of Fgf (Fig. 3(c)). Finally, Lfng KO simulation did not abolish the Hes7 oscillations, whereas the system oscillations were completely abolished in the Hes7 KO simulation (Fig. 3(d,e)). Hence, we concluded that the agreement of the logical model with experimental data is robust to the updating delays and to the initial conditions. Hence, the topology of this network and a set of relatively simple logical rules were able to reproduce PSM position dependent phenotypes of the Fgf and Notch mutants. 3.2. DDE model Although the previous logical model is the best approach to model genetic interactions where numerical parameter values are unknown, differential equation models might be more useful to gain insight in the numerical values of the biochemical reactions underlying the regulatory networks. Therefore, we next delineated a differential equation model, whose topology followed that of the previous logical model. The resulting DDE model included a delayed negative autoregulation of Hes7, and the Lfng mediated regulation of Notch: !   yrh h dh ln 2 nrn m Fgfr1 þ ð1  mÞ rn  h , ¼ th dt yn þ nrn yrh h þ hrT hh ! yrh h dl ln 2 nrn l , ¼ rn r tl yn þ nrn yh h þ hrT hh dt ! yrl l dn ln 2 n . (2) ¼ rl tn yl þ lrT ll dt The variables h; l; n represented the Hes7, Lfng and Notch1 activities, and hT h ¼ hðt  T h Þ; lT l ¼ lðt  T l Þ the delayed variables of Hes7 and Lfng. Those delays were estimated to be around

181

11.64–37 min (Hes7) and 10.5–60.2 min (Lfng) based on sequence data and around 37 min (Hes7 and Lfng) based on in situ hybridization data (Bessho et al., 2003) (Materials and methods). In the case of Notch1, the delay involving proteolysis and nuclear translocation should be short relative to the transcriptional regulation of Hes7 and Lfng. Hence, a good estimation of the Hes7 and Lfng transcription/translation delays is 37 min. However, a delay of 37 min resulted in small amplitude and damped oscillations (see below), so that we chose a delay of 45 min, which resulted in a reasonably good amplitude and was close to the theoretical estimations. To model the combined regulation of Hes7 by Fgf and Notch signaling, we considered an algebraic sum of the Fgfr1 and Notch1-dependent Hes7 synthesis to account for the data that the gene products (Notch1 and Fgfr1) independently induce the Hes7 promoter (Niwa et al., 2007). We chose Hill functions to describe the regulation of transcription and/or of synthesis of the molecules (including the regulation of Notch1 by the enzyme Lfng), because the sigmoidal shape of those functions made them very reminiscent of the sudden changes in transcription and/or synthesis observed in vivo. The normalized thresholds were yh ¼ kh ln ð2Þ=ðah th Þ; yl ¼ kh ln ð2Þ=ðal tl Þ; yn ¼ kh ln ð2Þ=ðan tn Þ (Materials and methods). From the parameters in this model, only the half-lives of Hes7 (th ¼ 22:3 min) and Lfng (tl ¼ 70 min) have been experimentally determined (Hirata et al., 2004; Shifley and Cole, 2008). Other parameter values have been estimated theoretically. The values of the synthesis rates ah ¼ al ¼ an ¼ 12 min1 , and thresholds values kh ¼ kl ¼ kn ¼ 4 (equivalent to the Michaelis-Menten constant) are in the order of previously proposed values (Lewis, 2003). Regarding the Hill coefficients, we chose a value 2 for Hes7 (rh ¼ 2), because Hes7 acts as a dimer. In the case of Lfng and Notch1, we chose values 1 (rl ¼ rn ¼ 1), because currently Lfng and Notch1 are not supposed to act as dimers. To model the Fgf signaling gradient, we could have defined an algebraic function. However, we preferred to use a model from the literature of the Fgf/RA module (Fig. 1(b)) to show that integration of different models is possible (Goldbeter et al., 2007). That model employed the reciprocal inhibition of the Fgf and RA gradients to achieve the sigmoidal shape of the Fgf signaling gradient (Goldbeter et al., 2007). We coupled this model to the model given by Eqs. (2) with minor notational changes (Materials and methods). Like for the logical model, the model was simulated in a single cell that moves between the posterior and anterior PSM. The time required by a new cell to reach the new somite was chosen to be 720 min based on arguments about the segmentation period, and PSM and new somite lengths (Materials and methods). It has been found that Fgf signaling acts in the posterior two-thirds of the PSM, i.e. 480 min in our simulations (Pourquie´, 2003). Then, at the beginning and at the end of the simulation (to480 min and t4480 min), one virtual cell perceived very high (posterior PSM) and very low (anterior PSM) level of Fgfr1 activation, respectively. The wild-type simulation showed that Hes7, Lfng and Notch1 oscillated in the posterior (high Fgfr1 values) and in the anterior PSM (low Fgfr1 values), which agreed with observations (Fig. 4(a)). To validate the model in regard to other experimental observations, we also carried out mutant simulations. In the case of the Fgfr1 cKO ðFgfr1 ¼ 0Þ, we observed that Hes7 quickly settled down to a stable fixed point (Fig. 4(b)). This agreed with the observation that Hes7 oscillations are completely lost in the Fgfr1 cKO embryos (Niwa et al., 2007). Next, we simulated the Rbpj loss-of-function mutant (dn=dt ¼ 0; nð0Þ ¼ 0). In this simulation, Hes7 oscillations were quickly damped in the domain of low Fgfr1 activation (Fig. 4(c)), which corresponds to the anterior PSM in agreement with experimental observations (Niwa et al., 2007).

ARTICLE IN PRESS 182

´lez, R. Kageyama / Journal of Theoretical Biology 259 (2009) 176–189 A. Gonza

Fig. 3. Strong connected components (SCCs) of the full state transition graph of wild-type and mutant situations for the logical parameters given in Eqs. (1). SCCs are useful to compress oscillatory regions of the state transition graph. SCCs were labeled as ‘‘cc-i’’ with the index i 2 f1; 2; . . .g and a green background color. States that did not belong to SCCs were labeled as ‘‘u-s’’ where s corresponded to the four-digits state of the system and had red (stable state) or gray background (other states). In the wild-type simulation (a), there were three SCCs (cc-0, cc-1 and cc-2) with high Fgf levels and one SCC (cc-3) with low Fgf levels meaning that the system oscillated under high and low Fgf levels (a). In the Fgfr1 cKO simulation (b), there was only some transient oscillations (cc-0) that corresponded to damping oscillations before reaching an stable steady state. In the Rbpj cKO simulation (c), the system oscillated only in the presence of high Fgfr1 levels (SCCs cc-0 and cc-1), i.e. in the posterior PSM. In the Lfng KO mutant simulation (d), the system oscillated in the presence of high (SCCs cc-0, cc-1 and cc-2) and low Fgfr1 levels (SCC cc-3) (d). Finally, in the Hes7 mutant simulation (e), the system did not oscillate as shown through the lack of SCCs. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The simulations of the Lfng mutation (dl=dt ¼ 0; lð0Þ ¼ 0) did not severely affect Hes7 oscillations (Fig. 4(d)). By contrast, simulation of a Hes7 loss-of-function mutation (dh=dt ¼ 0; hð0Þ ¼ 0) perturbed the oscillations of the Lfng expression as reported in the experimental literature (Bessho et al., 2001; Niwa et al., 2007).

In summary, our DDE model like the previous logical model was able to explain the PSM position dependent phenotypes of Hes7 in Fgf and Notch mutants. Taken together, the simulations of the DDE and the logical models suggested that the observed PSM dependent phenotypes can be explained though the negative

ARTICLE IN PRESS ´lez, R. Kageyama / Journal of Theoretical Biology 259 (2009) 176–189 A. Gonza

Fgfr1 cKO simulation

Wild-type simulation Fgfr1 RAR Hes7 Lfng Notch1

0.20

0.15

0.10

0.10

0.05

0.05

0.00

0.00 100

200

300 400 Time [min]

500

600

Fgfr1 RAR Hes7 Lfng Notch1

0.20

Levels

Levels

0.15

0

700

0

100

200

Rbpj cKO simulation

300 400 Time [min]

500

600

700

Lfng KO simulation Fgfr1 RAR Hes7 Lfng Notch1

0.20

Fgfr1 RAR Hes7 Lfng Notch1

0.20

0.15 Levels

0.15 Levels

183

0.10

0.05

0.10

0.05

0.00

0.00 0

100

200

300 400 Time [min]

500

600

700

0

100

200

300 400 Time [min]

500

600

700

Hes7 KO simulation Fgfr1 RAR Hes7 Lfng Notch1

0.20

Levels

0.15

0.10

0.05

0.00 0

100

200

300 400 Time [min]

500

600

700

Fig. 4. Wild-type (a) and mutant simulations (b–e) of the DDE model given by Eqs. (2). All simulations comprised a 720 min pre-simulation (to eliminate transients) followed by the plotted 720 min simulation. The x-axis showed the simulation time of 720 min, and indirectly the cell position in the PSM with the most posterior and most anterior position in the left and right side of the x-axis, respectively. When the simulation times were smaller than 480 min and larger than 480 min, the Fgfr1 activations were high and low, respectively. (a) In the wild-type simulation, Hes7 oscillated in the anterior and posterior regions. (b) In the Fgfr1 cKO mutant, Hes7 did not oscillate. (c) In the Rbpj cKO mutant (c), Hes7 only oscillated in the posterior PSM. (d and e) The Lfng mutation did not affect the oscillations in contrast to the Hes7 mutation that completely abolished the oscillations. Note that the curves of Fgfr1 and RAR were scaled down four times to fit them in the plots. Also note that in the Lfng KO simulation and Hes7 KO simulations (d and e), Notch1 and Lfng levels increased rapidly out of the plot. All simulations were initiated with hð0Þ ¼ 0:025; lð0Þ ¼ 0:6; nð0Þ ¼ 0:07; RAð0Þ ¼ 0; Mcð0Þ ¼ 7:45; Cð0Þ ¼ 26:62; Fð0Þ ¼ 5: In order to simulate the mutants, the equations were changed to: Fgfr1 ¼ 0 (b, Fgfr1 cKO); dn=dt ¼ 0 and nð0Þ ¼ 0 (c, Rbpj cKO); dl=dt ¼ 0 and lð0Þ ¼ 0 (d, Lfng KO) and dh=dt ¼ 0 and hð0Þ ¼ 0 (e, Hes7 KO).

ARTICLE IN PRESS 184

´lez, R. Kageyama / Journal of Theoretical Biology 259 (2009) 176–189 A. Gonza

feedback loop of Hes7 together with the double activation of Hes7 by the Fgf and Notch signalling pathways.

3.3. Robustness of the model to parameter changes Next, we analyzed how changes in parameters affected the oscillations. For this aim, we carried out numerical continuations where Fgfr1 was treated as a parameter (Figs. 5–7). We examined how the Hopf bifurcation was affected by the parameters of the DDE model. In the wild-type situation (m ¼ 0:8), the oscillations changed from sustained (Fgfr140:3) to damped (Fgfr1o0:3) (Fig. 5(a)). When the parameter m was decreased, the Hopf bifurcation moved anteriorly and eventually disappeared (Fig. 5(a)). This was due to the fact that for small m values, the Hes7 oscillations became independent of Fgf signaling, and oscillated all along the PSM independently of the Fgf signalling gradient. When the Hes7 delay T h was decreased, the region of damped oscillations moved towards the posterior PSM (right side of the x-axis) (Fig. 5(b)). By contrast, the Hopf bifurcation did not significantly change, when the Lfng delay T l was changed. The system oscillations were very sensitive to the Hes7 half-life in agreement with experimental data (Hirata et al., 2004) (Fig. 5(c)). By contrast, the Hopf bifurcation was relatively insensitive to large values of Lfng and Notch half-lives. This was in good agreement with measurements that showed a relative long Lfng protein halflife tl ¼ 70 min (Shifley and Cole, 2008). In the case of the synthesis rates, changes in the Lfng and Notch1 synthesis rates did not affect the oscillations (Fig. 5(d)). By contrast, a decrease of the synthesis rate, e.g. by using a transcription inhibitor, could result in damped oscillations. However, the sensitivity of the Hopf bifurcation to the Hes7 synthesis rate was much lower than to the Hes7 half-life as, for instance, decreasing by half the Hes7 synthesis rate did not completely perturb the oscillations. Changes of the Lfng and Notch1 thresholds (kl ; kn ) would not affect the oscillations (Fig. 5(e)). However, an increase of the Hes7 threshold kh could result in damped oscillations. Finally, changing the Hill coefficients of Lfng or Notch1 would not affect the oscillations so much (Fig. 5(f)). But even the slightest change in the Hill coefficient of Hes7 would completely abolish the oscillations. In our model, the cooperativity or Hill coefficient of Hes7 was the most sensitive parameter. We also examined how the oscillation period changed with parameter variations. According to the clock-and-wavefront model, the oscillation period sets the pace of segmentation and the length of the somites. On the other hand, it is unknown why the oscillation period increases in the anterior PSM. To gain insight into the parameters that might control the period, we have plotted the Hopf period as a function of different parameter values (Fig. 6). We found that the parameters that most affected the oscillation period were the Hes7 delay, half-life and very specially, the Hill coefficient (Figs. 6(b,c,f)). According to our model, larger values of those parameters would increase the oscillation period. Other changes in parameter values like larger Fgfr1 and Hes7 synthesis rates, or lower Hes7 threshold concentrations would also increase the period of the oscillations but to a much lesser extent (Figs. 6(a,d)). In the simulations of the DDE model, the amplitude and average concentrations seemed to largely change with different Fgfr1 levels (Fig. 4). Therefore, we also examined how changes in parameter values affected the amplitude of oscillations (Fig. 7). In these plots, the parameters that least affected the oscillations amplitude of all variables were the relative Fgfr1 contribution m and the Lfng delay T l (Fig. 7(b,d)). In a second group of parameters, the Hes7 solution profile was nearly constant, while the Lfng and Notch1 oscillation amplitudes and average concentration changed

to different extents. All the parameters of Lfng (except the Lfng delay T l ) and Notch1 belonged to this second group of parameters, i.e. half-lives (tl ; tn ), synthesis rates (al ; an ), thresholds (kl ; kn ) and Hill coefficients (rl ; rn ) (Fig. 7(f,g,i,j,l,m,o,p)). In a third group of parameters, the solution profiles of all variables Hes7, Lfng and Notch1 did significantly change. Fgfr1 and all Hes7 parameters belonged to this group (Fig. 7(a,c,e,h,k,n)). Thus, the Hes7 parameter changes affected very much the amplitude and average concentration of the oscillations of all variables, while the Lfng and Notch1 parameters did not affect so much the amplitude and average concentration of the Hes7 oscillations. Altogether, our results showed that the system is sensitive to variations of Hes7 parameter values and robust to changes in Lfng and Notch1 parameter values.

4. Discussion In this work, we wanted to formalize recent observations where PSM position-dependent phenotypes were reported for different mutants of Notch and Fgf signaling (Niwa et al., 2007). In the absence of Fgf signaling, Hes7 oscillations were lost after 2 h in the posterior and after 6 h also in the anterior PSM, whereas in the absence of Notch signaling, only anterior Hes7 oscillations were lost (Niwa et al., 2007). Our analysis provided explanation to these observations. In the absence of Fgf, Hes7 is inactive in a stable fixed point. After entering the PSM, cells come in contact with Fgf signaling, which destabilizes the Hes7 fixed point, and Hes7 starts to oscillate (Fig. 8). When the Fgf signaling levels become lower again, the oscillations undergo a Hopf bifurcation towards the stable fixed point. However, because the system is excited, these damped oscillations are maintained long enough for the cells to reach the forming somite through the Notch pathway in the anterior PSM. This explanation is based on the theoretical and computational analyses of a logical and a DDE models of the process in question. To reach these conclusions we delineated a model and challenged it with simulations of genetic perturbation experiments. Most of the observations could be reproduced. However, the results of some mutant simulations were valid only for specific developmental stages. In the Lfng KO situation, Hes7 expression oscillation was not affected (Fig. 4(d)). These results reproduces only what happens in vivo at embryonic stage E9.5 and later. Chen et al. (2005) reported that in E9.0 Lfng KO embryos, Hes7 expression oscillation is replaced by homogeneous expression in the PSM. By contrast, in E9.5 Lfng KO embryos, Niwa et al. (2007) reported that though Hes7 expression is affected, it still oscillates. Shifley et al. (2008) resolved this discrepancy by observing that in E8.5 Lfng mutant embryos, Hes7 shows homogeneous transcription, whereas in E10.5 Lfng KO embryos, Hes7 transcription oscillates. Therefore, our model should reproduce what is happening after the E9.5 stage. This underlined the importance to simulate mutant simulations to evaluate the model. Other mutant simulations defied our understanding of the system. Reports by some authors showed upregulation of Lfng expression in the Fgfr1 cKO (Niwa et al., 2007), while other authors reported downregulation of Lfng (Wahl et al., 2007). Our Fgfr1 cKO simulation reproduced the former observations. We tried to get the downregulation of Lfng by assuming that Lfng expression is also regulated by Fgf signaling (data not shown). Our simulations resulted in similar PSM position dependent phenotypes to Hes7, which has not been reported in Notch signaling mutants such as Dll1 (Chen et al., 2005), Dll3 (Chen et al., 2005), Psen1 (Koizumi et al., 2001) and after treatment with the gammasecretase inhibitor DAPT (Dale et al., 2006). Therefore, it remained

ARTICLE IN PRESS ´lez, R. Kageyama / Journal of Theoretical Biology 259 (2009) 176–189 A. Gonza

185

Fig. 5. Hopf bifurcation as a function of the DDE model parameters. For each equation parameter, we plotted the Hopf bifurcations in a two parameter space as a function of Fgfr1 (always in the x-axis) and the Fgfr1 contribution m (a), the delays (b), the half-lives (c), the synthesis rates (d), the threshold concentrations (e) and the Hill coefficients (f) in the y-axis. Parameter values on each side of the curve resulted in each either damped or sustained oscillations as stated within the plot. The Hopf bifurcations after continuation of Lfng and Notch1 parameters tended to be vertical meaning that the Hopf bifurcation was robust to changes in those Lfng and Notch1 parameters. By contrast, the curves of Hes7 parameters tended to be horizontal meaning that the Hopf bifurcation was sensitive to changes in Hes7 parameters.

ARTICLE IN PRESS 186

´lez, R. Kageyama / Journal of Theoretical Biology 259 (2009) 176–189 A. Gonza

Fig. 6. Hopf period as a function of the DDE model parameters. We plotted the Hopf period as a function of different values of each system parameter, that is the Fgfr1 level and Fgfr1 contribution m (a), the delays (b), the half-lives (c), the synthesis rates (d), the threshold concentrations (e) and the Hill coefficients (f). In all plots (except for Fgfr1 continuation in (a)), Fgfr1 was chosen to be value 0.8. The curves of the parameters of Lfng and Notch1 tended to be horizontal meaning that the Hopf period was robust to changes in these parameters. By contrast, the curves of the parameters of Hes7 tended to be more vertical meaning that the Hopf period was sensitive to changes in the Hes7 parameters.

ARTICLE IN PRESS ´lez, R. Kageyama / Journal of Theoretical Biology 259 (2009) 176–189 A. Gonza

187

Fig. 7. Amplitude of the oscillations as a function of the DDE model parameters. The two lines showed the maximal and minimal values of the oscillations of Hes7 (blue solid), Lfng (green dashed) and Notch1 (red dotted) as a function of the Fgfr1 levels (a), the Fgfr1 activation contribution m (b), the delays (c and d), the half-lives (e–g), the synthesis rates (h–j), the thresholds (k–m) and the Hill coefficients (n–p). In all plots (except (a)), Fgfr1 was chosen to be value 0.8. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

unclear why Lfng expression was downregulated in some Fgfr1 cKO embryos (Wahl et al., 2007). Our model allowed us to evaluate not so well understood aspects of the system like the qualitative, i.e. positive or negative, effect of Notch1 glycosylation by Lfng. We assumed a negative effect of this posttranscriptional modification based on the observation that in Lfng KO mouse embryos, the Notch1 intracellular domain is found throughout the PSM (Morimoto et al., 2005). Moreover, overexpression of Lfng in chick embryos downregulates Notch1 target genes (Dale et al., 2003). However, computer simulations indicated that the same phenotypes could be observed if we assume that Lfng sensitizes Notch1, and this

sensitized form has a shorter half-life (Cinquin, 2007). We simulated this hypothesis by assuming that Lfng activates Notch1 in the wild-type and Lfng mutant situation. Under this assumption, in the Lfng KO, Notch1 did not become activated, which does not agree with experimental observations. Therefore, we favored the assumption that Lfng negatively regulates the Notch receptor (Morimoto et al., 2005). Our modeling results were supported by two modeling formalisms, a logical and a DDE formalism. Combination of different modeling formalisms allowed us to take advantage of the strengths of each formalism. Also, this allowed us to confirm that our theoretical results did not depend on the specific

ARTICLE IN PRESS ´lez, R. Kageyama / Journal of Theoretical Biology 259 (2009) 176–189 A. Gonza

188

Fgf signaling

[Hes7]

Notch signaling Hopf bifurcation

Posterior

Anterior

Fig. 8. Summary of the conclusions of our computer simulations. Hes7 expression oscillations are induced by Fgf signaling in the posterior PSM. When the cell reaches the anterior end of the Fgf gradient, Hes7 oscillations become damped, i.e. there is a Hopf bifurcation. However, Hes7 expression continues to oscillate due to the initial excitation by Fgf signaling and the presence of Notch signaling.

bifurcation arises in the Hes1 network (Feng and Navaratna, 2007; Verdugo and Rand, 2008). However, our study is the first work to propose explicitly the induction of a Hopf bifurcation by Fgf signaling in the segmentation clock. In this work, we presented a network model of the mammalian segmentation clock under the control of the Fgf gradient in a single cell that moves towards the PSM. We showed with different modeling formalisms that striking phenotypes of Hes7 expression oscillations in Notch and Fgf mutants could be reproduced. Therefore, our model provides a significant contribution to our understanding of how the Fgf gradients control the segmentation clock.

Funding formalization but in the topology of the chosen network. As shown here, parameterization of the logical model was very straightforward to represent genetic interactions. However, differential equations gave us a finer representation of the system dynamics and allowed us to use powerful numerical continuation tools. Both formalisms agreed with most of the experimental observations as shown through wild-type and mutant simulations. Indeed, model predictions depend very much on the formalism. DDE models have been very successful in explaining the importance of transcription/translation delays (Lewis, 2003; Monk, 2003; Jensen et al., 2003) and the importance of a small enough Hes7 half-life (Hirata et al., 2004). By contrast, ordinary differential equation (ODE) models have been more successful in proposing explanations for the waves of gene expression that spread from the posterior to the anterior PSM (Tiedemann et al., 2007; Uriu et al., 2009). In our DDE model, we found that the period depends mainly on the Hes7 delay, half-life and Hill coefficient (Fig. 6). However, those parameters also affected very much the amplitude of the oscillations (Fig. 7). Also, the different levels of Fgfr1 affected very much the amplitudes and average concentrations of the model variables (compare high and low Fgfr1 levels in Fig. 4). To compensate the smaller amplitude of the oscillations under low Fgfr1 levels, we tested the possibility that other parameters changed with Fgfr1 or that the system was externally forced by another oscillator (data not shown). These different models could not correct the smaller amplitude in some cases or did not agree with mutant phenotypes in other cases. Therefore, a segmentation clock model is still to be found that is sensitive to the half-life, but where the period is controlled by some biochemical parameter without affecting the amplitude. Interestingly, it was found recently that by interlinking negative and positive feedback loops it is much easier to adjust the period without affecting the amplitude (Tsai et al., 2008). The idea of destabilization of Hes7 expression in response to the Fgf seems to be valid for the gene expression oscillation of Hes1 (member of Hes7 gene family) in fibroblasts. In fibroblasts, it was observed that the Hes1 expression starts to oscillate after serum shock (Hirata et al., 2002). Serum contains Fgf, and culture of fibroblasts with an Fgf inhibitor abolished these oscillations (Nakayama et al., 2008). This suggests that the same Fgf-dependent system destabilization that acts in the PSM and results in Hes7 oscillations, might also act in fibroblasts and result in Hes1 oscillations. Indeed, in fibroblasts, Hes1 expression oscillations were lost after treatment with the Fgf inhibitor SU5402 but were recovered after inhibitor washout (Nakayama et al., 2008). This suggests that after inhibitor washout, the expression oscillation of Hes1 was recovered again because of destabilization of the fixed point of Hes1. Recent theoretical models have studied the conditions under which the Hopf

Research in our laboratory was supported by the Genome Network Project; the Grants-in-Aid from the Ministry of Education, Culture and Sports, Science and Technology of Japan; and the Uehara Memorial Foundation. A.G. was supported by Postdoctoral Grant P06237 awarded by the Japan Society for the Promotion of Science.

Acknowledgments We thank Y. Niwa for helpful discussions and R. Szalai for advice with the PDDECONT software. Also, we acknowledge Y. Ma, D. Thieffry, A. Este´vez-Torres and C. Lamy for critical reading of this manuscript. References Audibert, A., Weil, D., Dautry, F., 2002. In vivo kinetics of mrna splicing and transport in mammalian cells. Mol. Cell. Biol. 22, 6706–6718. Aulehla, A., Wehrle, C., Brand-Saberi, B., Kemler, R., Gossler, A., Kanzler, B., Kanzler, B., Herrmann, B.G., 2003. Wnt3a plays a major role in the segmentation clock controlling somitogenesis. Dev. Cell 4, 395–406. Bessho, Y., Hirata, H., Masamizu, Y., Kageyama, R., 2003. Periodic repression by the bHLH factor Hes7 is an essential mechanism for the somite segmentation clock. Genes Dev. 17, 1451–1456. Bessho, Y., Sakata, R., Komatsu, S., Shiota, K., Yamada, S., Kageyama, R., 2001. Dynamic expression and essential functions of Hes7 in somite segmentation. Genes Dev. 15, 2642–2647. Chen, J., Kang, L., Zhang, N., 2005. Negative feedback loop formed by Lunatic fringe and Hes7 controls their oscillatory expression during somitogenesis. Genesis 43, 196–204. Cinquin, O., 2007. Repressor dimerization in the zebrafish somitogenesis clock. PLoS Comput. Biol. 3, e32. Cooke, J., Zeeman, E.C., 1976. A clock and wavefront model for control of the number of repeated structures during animal morphogenesis. J. Theor. Biol. 58, 455–476. Dale, J.K., Malapert, P., Chal, J., Vilhais-Neto, G., Maroto, M., Johnson, T., Jayasinghe, S., Trainor, P., Herrmann, B., Pourquie´, O., 2006. Oscillations of the snail genes in the presomitic mesoderm coordinate segmental patterning and morphogenesis in vertebrate somitogenesis. Dev. Cell 10, 355–366. Dale, J.K., Maroto, M., Dequeant, M.L., Malapert, P., McGrew, M., Pourquie´, O., 2003. Periodic Notch inhibition by Lunatic fringe underlies the chick segmentation clock. Nature 421, 275–278. del Barco Barrantes, I., Elia, A.J., Wu¨nsch, K., de Angelis, M.H.H., Mak, T.W., Rossant, J., Conlon, R.A., Gossler, A., de la Pompa, J.L., 1999. Interaction between Notch signalling and Lunatic fringe during somite boundary formation in the mouse. Curr. Biol. 9, 470–480. Deque´ant, M.-L., Glynn, E., Gaudenz, K., Wahl, M., Chen, J., Mushegian, A., Pourquie´, O., 2006. A complex oscillating network of signaling genes underlies the mouse segmentation clock. Science 314, 1595–1598. Dubrulle, J., Pourquie´, O., 2004. fgf8 mRNA decay establishes a gradient that couples axial elongation to patterning in the vertebrate embryo. Nature 427, 419–422. Engelborghs, K., Luzyanina, T., Roose, D., 2002. Numerical bifurcation analysis of delay differential equations using dde-biftool. ACM Trans. Math. Softw. 28, 1–21 ISSN 0098-3500. Ermentrout, B., 2002. Simulating, Analyzing, and Animating Dynamical Systems: A Guide Toi Xppaut for Researchers and Students. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA. ISBN 0898715067.

ARTICLE IN PRESS ´lez, R. Kageyama / Journal of Theoretical Biology 259 (2009) 176–189 A. Gonza

Feng, P., Navaratna, M., 2007. Modelling periodic oscillations during somitogenesis. Math. Biosci. Eng. 4, 661–673. Giudicelli, F., Ozbudak, E.M., Wright, G.J., Lewis, J., 2007. Setting the tempo in development: an investigation of the zebrafish somite clock mechanism. PLoS Biol. 5, e150. Goldbeter, A., Gonze, D., Pourquie, O., 2007. Sharp developmental thresholds defined through bistability by antagonistic gradients of retinoic acid and fgf signaling. Dev. Dyn. 236, 1495–1508. Goldbeter, A., Pourquie´, O., 2008. Modeling the segmentation clock as a network of coupled oscillations in the notch, wnt and fgf signaling pathways. J. Theor. Biol. 252, 574–585. Gonza´lez, A., Kageyama, R., 2007. Practical lessons from theoretical models about the somitogenesis. Gene Regul. Syst. Biol. 1, 35–42. Gonza´lez Gonza´lez, A., Naldi, A., Sanchez, L., Thieffry, D., Chaouiya, C., 2006. GINsim: a software suite for the qualitative modelling, simulation and analysis of regulatory networks. Biosystems 84, 91–100. Gridley, T., 2006. The long and short of it: somite formation in mice. Dev. Dyn. 235, 2330–2336. Hirata, H., Bessho, Y., Kokubu, H., Masamizu, Y., Yamada, S., Lewis, J., Kageyama, R., 2004. Instability of Hes7 protein is crucial for the somite segmentation clock. Nat. Genet. 36, 750–754. Hirata, H., Yoshiura, S., Ohtsuka, T., Bessho, Y., Harada, T., Yoshikawa, K., Kageyama, R., 2002. Oscillatory expression of the bHLH factor Hes1 regulated by a negative feedback loop. Science 298, 840–843. Jensen, M.H., Sneppen, K., Tiana, G., 2003. Sustained oscillations and time delays in gene expression of protein Hes1. FEBS Lett. 541, 176–177. Koizumi, K., Nakajima, M., Yuasa, S., Saga, Y., Sakai, T., Kuriyama, T., Shirasawa, T., Koseki, H., 2001. The role of presenilin 1 during somite segmentation. Development 128, 1391–1402. Lewis, J., 2003. Autoinhibition with transcriptional delay: a simple mechanism for the zebrafish somitogenesis oscillator. Curr. Biol. 13, 1398–1408. MacDonald, N., 1989. Biological Delay Systems: Linear Stability Analysis. Cambridge University Press, Cambridge. Monk, N.A., 2003. Oscillatory expression of Hes1, p53, and NF-kB driven by transcriptional time delays. Curr. Biol. 13, 1409–1413. Morimoto, M., Takahashi, Y., Endo, M., Saga, Y., 2005. The Mesp2 transcription factor establishes segmental borders by suppressing Notch activity. Nature 435, 354–359. Nakayama, K., Satoh, T., Igari, A., Kageyama, R., Nishida, E., 2008. FGF induces oscillations of Hes1 expression and Ras/ERK activation. Curr. Biol. 18, R332–R334. Niwa, Y., Masamizu, Y., Liu, T., Nakayama, R., Deng, C.-X., Kageyama, R., 2007. The initiation and propagation of Hes7 oscillation are cooperatively regulated by Fgf and Notch signaling in the somite segmentation clock. Dev. Cell 13, 298–304. O’Brien, T., Lis, J.T., 1993. Rapid changes in drosophila transcription after an instantaneous heat shock. Mol. Cell. Biol. 13, 3456–3463.

189

Pourquie´, O., 2003. The segmentation clock: converting embryonic time into spatial pattern. Science 301, 328–330. Rodriguez-Gonzalez, J.G., Santillan, M., Fowler, A.C., Mackey, M.C., 2007. The segmentation clock in mice: interaction between the wnt and notch signalling pathways. J. Theor. Biol. 248, 37–47. Santilla´n, M., Mackey, M.C., 2008. A proposed mechanism for the interaction of the segmentation clock and the determination front in somitogenesis. PLoS ONE 3, e1561. Serth, K., Schuster-Gossler, K., Cordes, R., Gossler, A., 2003. Transcriptional oscillation of lunatic fringe is essential for somitogenesis. Genes Dev. 17, 912–925. Shifley, E.T., Cole, S.E., 2008. Lunatic fringe protein processing by proprotein convertases may contribute to the short protein half-life in the segmentation clock. Biochim. Biophys. Acta. 1783, 2384–2390. Shifley, E.T., Vanhorn, K.M., Perez-Balaguer, A., Franklin, J.D., Weinstein, M., Cole, S.E., 2008. Oscillatory lunatic fringe activity is crucial for segmentation of the anterior but not posterior skeleton. Development 135, 899–908. Szalai, R., 2005. DDE-CONT: a continuation and bifurcation software for delaydifferential equations. Technical report, Budapest University of Technology and Economics, Hungary, 2005; available online from hhttp://seis.bris.ac.uk/ rs1909/pdde/i. Tam, P.P., 1981. The control of somitogenesis in mouse embryos. J. Embryol. Exp. Morphol. 65 (Suppl.), 103–128. Thieffry, D., Thomas, R., 1995. Dynamical behaviour of biological regulatory networks—II. Immunity control in bacteriophage lambda. Bull. Math. Biol. 57, 277–297. Thomas, R., D’Ari, R., 1990. Biological Feedback. CRC Press, Boca Raton, FL. Tiedemann, H.B., Schneltzer, E., Zeiser, S., Rubio-Aliaga, I., Wurst, W., Beckers, J., Przemeck, G.K., de Angelis, M.H., 2007. Cell-based simulation of dynamic expression patterns in the presomitic mesoderm. J. Theor. Biol. 248, 120–129. Tsai, T.Y.-C., Choi, Y.S., Ma, W., Pomerening, J.R., Tang, C., Ferrell, J.E., 2008. Robust, tunable biological oscillations from interlinked positive and negative feedback loops. Science 321, 126–129. Ujvari, A., Aron, R., Eisenhaure, T., Cheng, E., Parag, H.A., Smicun, Y., Halaban, R., Hebert, D.N., 2001. Translation rate of human tyrosinase determines its nlinked glycosylation level. J. Biol. Chem. 276, 5924–5931. Uriu, K., Morishita, Y., Iwasa, Y., 2009. Traveling wave formation in vertebrate segmentation. J. Theor. Biol. 257, 385–396. Verdugo, A., Rand, R., 2008. Hopf bifurcation in dde model of gene expression. Commun. Nonlinear Sci. Numer. Simulation 13, 235–242. Wahl, M.B., Deng, C., Lewandoski, M., Pourquie´, O., 2007. FGF signaling acts upstream of the NOTCH and WNT signaling pathways to control segmentation clock oscillations in mouse somitogenesis. Development 134, 4033–4041. Zhu, L., Hu, C., Li, J., Xue, P., He, X., Ge, C., Qin, W., Yao, G., Gu, J., 2007. Real-time imaging nuclear translocation of akt1 in hcc cells. Biochem. Biophys. Res. Commun. 356, 1038–1043.