Short-lived Her proteins drive robust synchronized oscillations in the zebrafish segmentation clock

Oscillatory systems pervade biology on every scale. Predator-prey population dynamics, daily rhythms of organisms and developmental pattern formation all show oscillatory behavior over periods of time. The anterior-posterior body axis of vertebrates is patterned as a series of spatially repetitive somite segments. As cells emerge from the presomitic mesoderm (PSM), located at the growing posterior end of the embryo, they are grouped together and organized into segments known as somites that will later differentiate into the vertebral column and trunk, and limb muscle tissues (Pourquié, 2011).

The periodicity of somite segmentation is controlled by a gene-expression oscillator, called the segmentation clock, which ticks in the cells of the unsegmented mesoderm. The oscillation period of the segmentation clock dictates the period of somite segmentation. Breakdown of oscillations disrupts somite segmentation and results in vertebral defects (Pourquié, 2011). Somite segmentation, i.e. somitogenesis, occurs in all vertebrates and has been experimentally studied in many species (Gomez et al., 2008; Eckalbar et al., 2012).

Hairy/enhancer-of-split-related genes (called Hes/her genes) are the only conserved oscillating genes in all studied species (Krol et al., 2011; Eckalbar et al., 2012). In zebrafish, the genes thought to be involved in somitogenesis that show oscillatory expression in the posterior PSM include her1, her7 and deltaC (Ozbudak and Pourquié, 2008) (Fig. 1). her1 and her7 genes encode members of the basic helix-loop-helix (bHLH) family of transcription factors; they function as transcriptional repressors only following dimerization. deltaC encodes a ligand for the Notch receptor. Hes6 protein (also called Her13.2) belongs to the same bHLH family of repressors as Her1 and Her7, but its expression does not oscillate, unlike Her1 and Her7 (Kawamura et al., 2005). Her1, Her7 and Hes6 proteins form homo- and heterodimers at different levels (Schröter et al., 2012; Trofka et al., 2012; Hanisch et al., 2013), and these dimers repress transcriptions of her1, her7 and deltaC (Giudicelli et al., 2007) forming a negative-feedback loop that has the potential to create oscillatory gene expression. It is the oscillatory expression patterns of these transcriptional repressors within the PSM that are proposed to be the mechanism underlying the segmentation clock (Lewis, 2003). Each oscillatory cycle defines a new somite and individual cellular positioning within that somite depends on the corresponding phase of the cycle during which each cell was generated. Delta-Notch signaling enhances the transcription of her1 and her7 and ensures their oscillations are synchronized across neighboring cells (Jiang et al., 2000; Horikawa et al., 2006; Mara et al., 2007; Riedel-Kruse et al., 2007; Ozbudak and Lewis, 2008; Delaune et al., 2012) (Fig. 1).

Mathematical models have been employed frequently to obtain a mechanistic understanding of this fascinating developmental clock and shed light on its counterparts in different organisms. Ozbudak and Lewis (Ozbudak and Lewis, 2008) modeled the system using only her1, her7 and deltaC genes in a two-cell system at the posterior PSM. Cinquin (Cinquin, 2007) modeled the same system using her1, her7, hes6 and deltaC genes in a one-dimensional chain of cells. Both of these models were built on experimental data that have been extensively updated recently. The system has been recently modeled by Schröter et al. (Schröter et al., 2012) and Hanisch et al. (Hanisch et al., 2013). Although the model published by Schröter et al. matched some current experimental findings, it is a simplified deterministic model, concentrating only on a single cell and approximating dimerization of proteins rather than explicitly including dimer species. Owing to its simplicity, this model could not address the phenotypes observed with loss of Notch signaling (such as loss of synchrony among oscillations of neighboring cells, reduction in oscillation amplitude and increase in oscillation period). Furthermore, this model could not predict how synchronization is affected in the mutants of her family genes. Hanisch et al. modeled the system at the multicellular level, but this model did not address the mutant phenotypes of all her family genes (such as the hes6^-/- and her7^-/-;hes6^-/- mutants) and the change in the period of segmentation when Notch signaling was impaired.

In the present study, we extended these earlier models by building a multicellular stochastic computational model for the oscillations taking place at the posterior PSM. Our simulations demonstrated that a modified version of a mathematical model based on a simple autoregulatory feedback loop (Lewis, 2003) can accommodate the experimental observations (Jiang et al., 2000; Horikawa et al., 2006; Mara et al., 2007; Riedel-Kruse et al., 2007; Ozbudak and Lewis, 2008; Herrgen et al., 2010; Delaune et al., 2012; Schröter et al., 2012; Trofka et al., 2012; Hanisch et al., 2013). The autoinhibitory negative-feedback loops composed of two redundantly acting repressor dimers are sufficient to generate synchronized oscillations in wild-type embryos and we found parameter ranges that reproduce how the dynamics (i.e. period, amplitude and synchronization) of the segmentation oscillator are modified in different genetic backgrounds.

Previously, it has been proposed that the oscillating Her/Hes proteins should be short-lived for sustained oscillations to occur (Lewis, 2003). The half-life of mouse Hes7 protein was measured as 22 minutes in cell culture. Mutations that increase half-life of Hes7 disrupted oscillations and resulted in vertebral defects in mice (Hirata et al., 2004). Using our model, we randomly varied the biochemical parameters in biologically reasonable ranges and identified parameter sets that result in sustained synchronized oscillations in wild-type embryos and respective period, amplitude or synchronization changes in various mutant conditions. In all the identified parameter sets, half-life of Her proteins were shorter than 6.2 minutes, a rate that is roughly one quarter of the measured degradation rate in mice. We tested the prediction of our model by measuring the half-life of Her7 protein in a transgenic animal that can express an HA-tagged version of Her7 protein in an inducible manner. We found the half-life of HA-Her7 to be 3.5 minutes, matching the prediction of our model.

We have developed a time-delayed differential equation (DDE) model consisting of 14 equations with 44 parameters. Each equation describes the rate of change of an mRNA, monomer protein or dimer protein.

The variables: mh_i where i ∊ {1, 7, 6} and md represent the number of mRNA molecules of her1, her7, hes6 and deltaC, respectively; ph_i where i ∊ {1, 7, 6} and pd represent the number of protein molecules of Her1, Her7, Hes6 and DeltaC, respectively; and ph_i,j where i ≤ j and i,j ∊ {1, 7, 6} represent the number of molecules of Her1-Her1, Her1-Her7, Her1-Hes6, Her7-Her7, Her7-Hes6 and Hes6-Hes6 dimers. Each term in the model represents a rate of change for a biological reaction. Protein synthesis and degradation rates are described by psh_i, pdh_i where i ∊ {1, 7, 6}, psd and pdd for Her1, Her7, Hes6 and DeltaC, respectively. mRNA synthesis and degradation rates are described by msh_i, mdh_i where i ∊ {1, 7, 6}, msd and mdd for her1, her7, hes6 and deltaC, respectively. Dimer association, dissociation and degradation rates for Her1-Her1, Her1-Her7, Her1-Hes6, Her7-Her7, Her7-Hes6 and Hes6-Hes6 are represented by dah_i,j, ddh_i,j and pdh_i,j where i ≤ j and i,j ∊ {1, 7, 6}, respectively. DNA-binding dissociation constraints are critph_1,1, critph_6,7 and critpd, for Her1-Her1, Her7-Hes6 and Notch (NICD), respectively. Transcriptional and translational delays are represented by nmh_i where i ∊ {1, 7}, nph_i where i ∊ {1, 7, 6}, nmd and npd for her1, her7, hes6 and deltaC, respectively. In the model equations below, we represent the k^th cell as c_k and time as t.

Detailed descriptions of these equations and their derivations are provided in supplementary material Appendix S1.

The stochastic equations have been written as a direct translation of the deterministic model described above (supplementary material Appendix S1). All reactions were executed following an exponential law for which the reaction rate was the same as its deterministic counterpart. Similar to the deterministic model, protein-DNA binding and unbinding events (very fast compared with other events) have not been explicitly modeled. This is because reactions occurring on slower time scales (such as synthesis and degradation of mRNAs and proteins) smoothen the noise stemming from the reactions with fast dynamics (such as protein-DNA interactions) (Paulsson, 2004).

The deterministic simulation in our study is used to find biologically realistic parameter sets. The deterministic model is solved numerically using Euler’s method. Euler’s method increments the time in the chosen step size (e.g. 0.01 minutes), and updates mRNA and protein levels at each iteration using the rate of changes provided by the model. The stochastic simulations in our study are performed using the next reaction method (NRM), which discretely computes concentration levels based on probabilistic calculations (Anderson, 2007). Probabilistically determined propensities and reaction times are used to decide which reaction fires at each iteration. Reactions with higher propensities are more likely to fire. A delayed reaction queue is incorporated into the standard NRM algorithm to accommodate time delays (Anderson, 2007). Each iteration in NRM is computed as follows:

1. Update the propensity values related to the most recently fired reaction for each cell.

2. Calculate the time gap (the size of the next time step) using propensities.

3. Increment the time step and the relevant molecular counts.

4. If a delayed reaction is initiated, add it to the appropriate list. Otherwise fire immediate reactions and delayed reactions that are finished.

5. Repeat until simulation time expires.

Upper limits for the mRNA decay rates and transcriptional time-delays were measured previously (Giudicelli et al., 2007). However, most of the reaction rates in the segmentation clock network have not been experimentally measured precisely owing to technical difficulties. We have employed a parameter search to identify parameter values that were capable of reproducing experimental observations.

Initial parameter sets were composed of randomly generated parameter values (within biologically relevant ranges) (supplementary material Table S1) and deterministic simulations were then run with the resulting parameters. Initial parameter ranges using deterministic simulations have been constrained by comparison with experimental period observations for wild type, and her1^-/-, her7^-/-, hes6^-/-, her7^-/-;hes6^-/- and notch1a^-/- mutant zebrafish embryos (supplementary material Table S1). Approximately 5,000,000 parameter sets were run, with 200,000 passing the initial criteria on period restrictions. These 200,000 parameter sets were run with the NRM stochastic algorithm. The same restrictive period values were used to find the parameter sets that also worked for the two-cell stochastic simulations, resulting in a smaller number of parameter sets. The parameter values were further pruned by comparison to the experimental data on synchronization and amplitude of oscillations using two-cell stochastic simulations. Finally, the parameter sets that passed all of the previous steps (∼200) were used to run our stochastic simulation code for a sixteen-cell tissue (4×4). The simulations are compared with the current experimental data on period, synchronization and amplitude of oscillations in different genetic backgrounds. Only 41 parameter sets have passed these final criteria.

The last peak and trough of the her1 mRNA oscillations were used to calculate the period and amplitude of the oscillations. Period was calculated as the time difference between the last two peaks, and amplitude was calculated as the expression level change between the last peak and trough.

Calculation of oscillation features in the stochastic simulations was performed slightly differently owing to the noise in the resulting data. The first step in the calculation was smoothing the data using a 40 time step moving average method. Keeping the smoothing interval small did not affect the accuracy of the results. Because the peaks and troughs vary across the data, all the local maxima and minima were used as peak and trough values. A value more than five values to its left and five to its right was considered a local maximum, and a value less than five values to its left and five to its right was considered a local minimum. For every peak-trough pair, the value of the period and amplitude was calculated. To obtain an overall value of the period and amplitude for a run, these values were then averaged, and then averaged again over the number of runs when doing multiple runs.

In the stochastic simulations, the measure of synchronization was calculated using Pearson’s correlation coefficient. Within one run, the score was the average of the correlation coefficients between each cell and the first cell. For multiple runs, the average of the individual scores for each run was used.

We have used Student’s t-test to calculate P values and to determine the statistical significance of our results throughout the paper.

The codes for the study have been implemented in C++ and Python (available upon request). C++ was used because of its speed and Python because of its easy-to-use plotting libraries. The current version of our code can finish 1200 minutes of deterministic and stochastic simulations for a 16-cell tissue in less than a minute and four hours, respectively (on a node containing two quad-core Xeons at 2.3 GHz, 8 GB of memory). Parallel versions of the code have been written in order to run time-intensive stochastic simulations in parallel in a scientific cluster.

Heterozygous Tg(hsp70l:HA-her7) (Giudicelli et al., 2007) zebrafish embryos of ∼12-14 somite stages were dechorionated manually and heat shocked in a water bath for 1 hour at 37°C. Animals were then placed back in fish water at 24°C to recover for 5 minutes and then transferred to embryo media containing 100 μg ml^-1 cycloheximide for 15 minutes of incubation. For each measuring time point, 15 embryos were collected every 2.5 minutes and deyolked in 1 ml ice cold 1×Ringer’s solution without calcium using a glass Pasteur pipette. Embryos were then harvested by centrifugation at 350 g for 30 seconds at room temperature and supernatant was carefully taken off. The resulting cellular fractions were suspended in 48 μl protein lysis buffer [50 mM Tris-HCl, pH 7.4, 100 mM NaCl, 5 mM MgCl₂, 0.5 mM EDTA, 0.1% Triton X-100, supplemented with 10×EDTA-free protease inhibitor cocktail (Roche), 10 mg ml^-1 pepstatin, 10 mM sodium pyrophosphate pH 7.5, 10 mM β-glycerophosphate pH 7.5, 10 mM NaF] containing 250 U ml^-1 nuclease (Thermo Fisher) and subsequently lysed by adding 12 μl 5×SDS Laemmli protein sample buffer. Samples were immediately frozen in liquid nitrogen for storage.

For western blot analysis, the total protein extracts were heated at 85°C for 5 minutes and loaded at volumes of 37.5 μl (corresponding approximately to nine embryos) on 4-15% gradient gels and blotted overnight to PVDF membranes. Immunodetection of HA-Her7 proteins was carried out using mouse anti-HA (Roche 12CA5) and HRP-conjugated anti-mouse secondary antibodies. Relative amounts of detected HA-Her7 proteins were quantified after imaging using ImageJ (NIH) software.

Biological systems are often analyzed using deterministic simulations in which the input variables entirely determine the output. However, biological reactions are inherently noisy (stochastic) (Ozbudak et al., 2002). Stochastic simulations reflect this randomness by allowing random fluctuations to influence the output of biological processes. Thus, biological systems should be analyzed with stochastic models to capture a faithful representation of the systems. Furthermore, gene expression noise poses a great challenge to maintain the synchrony among oscillating neighboring cells in the segmentation clock (Riedel-Kruse et al., 2007; Ozbudak and Lewis, 2008; Delaune et al., 2012). Therefore, to obtain a more accurate description of this system and to determine how the synchronization of oscillations is affected when different genes in the segmentation clock network were mutated, we have utilized stochastic simulations.

Previous studies demonstrate that double knockdown of two oscillating genes, her1 and her7, prevents somite segmentation (Henry et al., 2002; Oates and Ho, 2002; Gajewski et al., 2003). The period of segmentation does not change in the single her1^-/- or her7^-/- mutants (Schröter et al., 2012; Hanisch et al., 2013). hes6^-/- and her7^-/-;hes6^-/- mutants do not have any segmentation defect (Sieger et al., 2006; Schröter et al., 2012; Trofka et al., 2012), but the period of the oscillations is increased by ∼6% (Schröter and Oates, 2010; Schröter et al., 2012). Double knockdown of her1 and hes6 also prevents somite segmentation. Mutations of genes involved in Notch signaling increases the period of segmentation by ∼7-23% over wild type (Herrgen et al., 2010). Our stochastic model demonstrates that the current zebrafish segmentation clock network is able to explain all quantitative experimental observations in wild type and her1^-/-, her7^-/-, hes6^-/-, her7^-/-;hes6^-/- and notch1a^-/- mutant embryos. We present synchronization (Fig. 2), period (Fig. 3) and amplitude (Fig. 4) values obtained for the her1 mRNA levels by our model. These values are similar for her7 mRNA levels.

Our model predicts a period of ∼31-32 minutes for her1 mRNA oscillations in wild-type embryos, a similar period in her1^-/- and her7^-/- (P=1 and 0.1, respectively), ∼5.7% increase in hes6^-/- (P=5×10^-8), 5.8% increase in her7^-/-;hes6^-/- (P=5.8×10^-9) and 8.5% increase in notch1a^-/- (P=6.5×10^-8) mutant embryos (Figs 3, 5). her1 mRNA oscillations are lost in her1^-/-;her7^-/- and her1^-/-;hes6^-/- double mutants. This is an obvious consequence of our model in which only Her1-Her1 and Her7-Hes6 dimers repress transcription of her1, her7 and deltaC. These are all in agreement with current experimental findings.

The segmentation clock is synchronized across neighboring cells via the Delta-Notch signal transduction pathway (Jiang et al., 2000; Horikawa et al., 2006; Mara et al., 2007; Riedel-Kruse et al., 2007; Ozbudak and Lewis, 2008; Delaune et al., 2012); synchronized gene expression is lost among neighboring cells when Notch signaling is impaired, which has been reproduced by our model (P=4×10^-31) (Figs 2, 5; supplementary material Movies 1-6).

her7^-/- mutant embryos exhibit segmentation defects that are confined to the posterior axis (Choorapoikayil et al., 2012; Schröter et al., 2012; Hanisch et al., 2013). Our model shows that synchronized gene expression among cells has been impaired in these embryos in comparison with wild type (P=2.8×10^-19) (Figs 2, 5). Somites segment properly in her1^-/-, hes6^-/- and her7^-/-;hes6^-/- mutants (Sieger et al., 2006; Schröter et al., 2012; Trofka et al., 2012; Hanisch et al., 2013). Our simulations show that synchronized gene expression among cells in her1^-/-, hes6^-/- and her7^-/-;hes6^-/- mutants is preserved, similar to wild-type embryos (P=0.5, 0.9 and 0.1, respectively) (Figs 2, 5).

A critical amplitude of gene expression oscillations is important for proper segment formation. It has been shown that, when Notch signaling is blocked with the inhibitor DAPT, the amplitude of her1 mRNA expression is decreased by ∼20% (Ozbudak and Lewis, 2008). In close agreement, our model provides an average decrease of ∼17% (P=6.3×10^-5) (Figs 4, 5). The amplitude of her1 mRNA expression is similar in her1^-/- mutant embryos to that in the wild-type embryos, which is reproduced in our model (P=0.7). Our model provides an explanation for the observation her7^-/- mutant embryos show defects in somite formation, suggesting that this might be due to reduction in amplitude and partial loss of synchronization among neighboring cells. Our simulations show that the amplitude of her1 mRNA expression decreases by ∼21% in her7^-/- mutant embryos in comparison with wild-type embryos (P=4.3×10^-7) (Figs 4, 5). Forty-one parameter sets satisfied all the experimentally observed changes in the period, amplitude or synchronization of oscillations in all mutant conditions. In 29 of these parameter sets, the amplitude of oscillations in single cells is reduced by >15% in her7^-/- mutants, in 31 of them synchronization is lost by >15% in her7^-/- mutants, and in 26 of them amplitude is reduced and synchronization is decreased by >15% in her7^-/- mutants in comparison with wild-type embryos. Experimental studies suggest that the amplitude of her1 mRNA expression is similar in hes6^-/- and her7^-/-;hes6^-/- mutant embryos to that in wild-type embryos (Schröter et al., 2012); our model reproduces this result with only an insignificant increase in amplitude: ∼2% in hes6^-/- and ∼3% in her7^-/-;hes6^-/- mutants (P=0.6 and 0.4, respectively) (Fig. 5).

The existence of sustained oscillations depends on the half-lives of oscillating mRNA and proteins (Lewis, 2003). An upper limit of 6-8 minutes was previously calculated for the half-lives of her1, her7 and deltaC mRNAs (Giudicelli et al., 2007). In our parameter searches, we obtained half-life parameters of between 2 and 6 minutes for these three mRNAs that satisfy the dynamic features of oscillations in all assessed genetic conditions. We calculated transcriptional time delays of 9-10 minutes for her1, her7 and deltaC in our recent study (Hanisch et al., 2013). Here, we obtained parameters corresponding to a range of 6-12 minutes for the transcriptional time delays for these three genes to satisfy all the experimental conditions. Overall, our simulation results agree with the experimentally measured or estimated reaction rates. However, future studies that more precisely measure these biochemical rates as well as other so far unmeasured rates are needed.

To measure the third critical parameter, the half-life of Her proteins, we utilized a transgenic animal that can inducibly express HA-Her7 upon heat-shock treatment of embryos (Giudicelli et al., 2007). We administered an hour of heat shock, recovered embryos for 5 minutes at 24°C and incubated them with 100 μg/ml cycloheximide for 15 minutes before collecting embryos every 2.5 minutes. We then performed western blots using anti-HA antibodies to quantify Her7 abundance and used an exponential decay curve to determine the degradation rate of HA-Her7 protein. We measured the half-life of HA-Her7 as 3.5 minutes at 24°C (Fig. 6B). This rate is consistent with our simulation parameters, where we find the half-life of Her proteins between 1.8 and 6.2 minutes (average=3 minutes, standard deviation=1.1 minutes) (Fig. 6A).

Somitogenesis is one of the best examples demonstrating the importance of developmental timing, and understanding this process in zebrafish lays the foundation for understanding it in more complex organisms. Experimental data have been obtained more quantitatively in recent years. We had previously measured the transcriptional time delays (Giudicelli et al., 2007; Hanisch et al., 2013) and obtained upper limits for the half-lives of her1, her7 and deltaC mRNAs (Giudicelli et al., 2007). We further demonstrated that the amplitude of her1 expression was reduced by ∼20% when Notch signaling was blocked by treatment of the inhibitor DAPT (Ozbudak and Lewis, 2008). Recent measurements of changes in the segmentation period in her1^-/-, her7^-/- (Schröter et al., 2012; Hanisch et al., 2013), hes6^-/- and her7^-/-;hes6^-/- (Schröter and Oates, 2010; Schröter et al., 2012) mutants provided quantitative data that could be utilized to build more realistic computational models.

In this study, we have modeled the zebrafish segmentation clock by focusing on the posterior PSM and using current knowledge of the segmentation clock network. Our stochastic multicellular model demonstrated that: (1) synchronized oscillations can be obtained in wild-type conditions; (2) the synchrony will be lost, the amplitude of oscillation will be reduced and the period of oscillation will be lengthened in notch1a^-/- mutants; (3) the period of oscillation, the amplitude of oscillation and the synchronization of oscillations are not affected in her1^-/- mutants; (4) oscillations are more sensitive to loss of her7 as the amplitude and the synchronization of oscillations are reduced in her7^-/- mutants; (5) the period of oscillations is not affected in her7^-/- mutants; (6) the period of oscillations is lengthened in hes6^-/- mutants; (7) the amplitude and the synchronization of oscillations are not affected in hes6^-/- mutants; (8) the period change in the her7^-/-;hes6^-/- mutant is similar to that of the hes6^-/- mutant; (9) sensitization of the her7^-/- mutant can be rescued by obtaining double her7^-/-;hes6^-/- mutants as the amplitude and the synchronization of oscillations are not affected in her7^-/-;hes6^-/- mutants; and (10) oscillations are abolished in her1^-/-;her7^-/- or her1^-/-;hes6^-/- double mutants.

We and others have shown that some of the posterior somites segment but other somites fail to segment in her7^-/- mutants (Choorapoikayil et al., 2012; Schröter et al., 2012; Hanisch et al., 2013). Our multicellular stochastic model captures this aspect accurately as it shows that the segmentation clock is sensitized against ‘noise’ in gene expression (Ozbudak et al., 2002) rather than damped out as previously claimed (Schröter et al., 2012). The synchronous oscillations are compromised; the synchrony levels fall below a threshold level during some of the oscillation cycles in a given embryo. But why does this happen? Figure 7A shows that the average DeltaC protein levels in her7^-/- mutants are much higher than those in wild type and in other her/hes mutant backgrounds. This is due to deltaC being de-repressed during the ‘off’ phase of sensitized oscillations. This deficiency results in more-or-less continuous (small peak-to-trough) Delta-Notch signaling between the neighboring cells. In other words, cells fail to deliver information regarding their oscillation phase to their neighboring cells. Furthermore, in her7^-/- mutants, more Hes6 proteins are free to dimerize with Her1 protein. This essentially creates an excess of Her1-Hes6 heterodimers that cannot bind to DNA targets (Schröter et al., 2012; Trofka et al., 2012). As a result, the total level of dimers that can repress transcription is reduced in comparison with wild type and other her/hes mutants (Fig. 7B). The amplitudes are in general lower because of the titration affect, as has been pointed out before (Schröter et al., 2012).

Binding of Hes6 increases the ‘effective’ degradation rates of Her1 and Her7 monomers. In other words, there are fewer Her1 or Her7 monomers available in the cell. Although the half-lives of Her proteins might not change in hes6^-/- and her7^-/-;hes6^-/- mutants, the outcome is equivalent to an increase in their half-lives, i.e. effectively there are more free Her1 (and Her7) proteins. Increasing the half-lives of oscillating proteins increases the period of oscillations (Lewis, 2003). The impact of Hes6 in reducing the levels of free Her monomers is strong because it is expressed continuously; hence, its average level is expected to be higher than those of Her1 and Her7. Our model agrees with previous assertions (Schröter et al., 2012) and shows that the effective levels of Her1 and Her7 are decreased in Hes6 and Her7;Hes6 mutants compared with other mutants (Fig. 7C; supplementary material Fig. S1).

Lewis (Lewis, 2003) derived a closed-form equation describing the period of oscillations in a single-cell oscillator as being mainly determined by the time delays and half-lives. How the period of oscillations in the multicellular system changes with respect to the Notch-Delta synthesis rate remained to be determined. Here, we obtained a relationship between the period of oscillations and the rate of synthesis of Notch-Delta proteins by numerically fitting the results obtained in our model. Our results show that the period of oscillations versus the DeltaC translation rate fits to an exponentially decreasing function (Fig. 7D).

Briefly, our analysis goes beyond that of Hanisch et al. (Hanisch et al., 2013) in demonstrating how the period of segmentation changes when Notch signaling is impaired and how the period, synchronization and amplitude of oscillations are affected in hes6^-/- and her7^-/-;hes6^-/- mutants. We also go beyond the analysis of Schröter et al. (Schröter et al., 2012) in showing that synchronized robust oscillations can be generated in wild-type embryos by coupling oscillations of redundant repressor dimers in neighboring cells via Notch signaling. We also show how the period and amplitude of oscillations are affected in notch1a^-/- mutant embryos, and how the synchronization of oscillations is affected in single and double mutants of her family genes.

The predictive capacity of a mathematical model is assessed by how well it can explain all the experimental data. In different biological systems, initial models developed to explain early experimental observations are often abandoned by the accumulation of additional data. Likewise, the rapid increase in the amount of experimental data on somite segmentation places constraints on the success of mathematical models. By matching to all quantitative data obtained in various genetic backgrounds so far, our model strongly supports the pacemaker model for the zebrafish somite segmentation clock, which is based on autoregulatory negative-feedback loops of short-lived Her proteins.

In addition to matching recently obtained quantitative data, our model also makes further predictions for future experiments. For example, the amplitude changes in individual or double her mutants and the compromised synchrony in her7^-/- mutants can be measured by using transgenic animals reporting oscillations in real time (Delaune et al., 2012) to test the predictions of our model.

There are important experimental questions that have not yet been quantitatively assessed and which could limit the predictive capacity of any computational model. First, the two DNA-binding Her dimer species (Her1-Her1 and Her7-Hes6) might not repress all the oscillating genes, particularly her1, her7 and deltaC, with equal strength. DNA binding does not guarantee efficient repression, which might also depend on cooperativity among binding sites for repressors and competition among different DNA-binding transcription factors. Second, our model also excludes a number of other her family genes such as her4, her12 and her15. Knockdown of these genes does not impact the segmentation process. However, they seem to contribute to stripy expression of oscillating genes in the anterior PSM (Shankaran et al., 2007). Hence, they might play auxiliary roles in segmentation. Generation of mutants for these genes and crossing them with her1 or her7 mutants would assess whether any of these genes contribute to the functional repressor pools in the system. Third, Delta proteins might dimerize and each homo- or heterodimer might have different activating or inhibitory roles on Notch signaling (Wright et al., 2011). This might explain the differences in the change of period among different Notch pathway mutants (Herrgen et al., 2010). Future quantitative data can be incorporated into our model to elucidate the differential changes of segmentation period in different mutants of genes influencing Notch signaling. Owing to these shortcomings, diverse parameter sets might be able to produce the same qualitative observations. However, if more quantitative data are produced, our model can be updated to narrow down the parameter ranges. As quantitative data are completely lacking for the processes ongoing in the anterior PSM, we chose to focus on modeling the segmentation clock in the posterior PSM. We plan to expand our model to accommodate processes taking place in the anterior PSM as quantitative data are produced.

One of the key parameters that strongly influences the existence and period of oscillations is the half-life of oscillating Her proteins. This value had so far been missing in zebrafish, but was measured to be ∼22 minutes for mouse Hes7 protein (Hirata et al., 2004). As the clock speed is much faster in zebrafish, even a half-life of 22 minutes will not allow oscillations to occur in zebrafish. Via exhaustive parameter searches, our model predicts an average half-life of ∼3 minutes for Her proteins; occasionally we obtained half-lives of up to 6.2 minutes for Her7 protein that were able to satisfy the quantitative features of somite segmentation in all mutant backgrounds (Fig. 6A).

We utilized a transgenic animal in which HA-tagged Her7 (HA-Her7) can be inducibly expressed in zebrafish embryos (Giudicelli et al., 2007). Previously, HA-Her7 was demonstrated to decay at a near-uniform rate in all tissues in the embryo, suggesting that the half-life of Her7 could be roughly equal in all tissues. We induced expression of HA-Her7 in the embryos and treated embryos with cycloheximide (to arrest translation) and measured the half-life of HA-Her7 as 3.5 minutes (Fig. 6B). We noticed that the initial rapid decay in Her7 protein levels saturates at later time points. There are two potential causes of this observation. (1) If there is cooperative stability (Buchler et al., 2005), i.e. dimers decay slower than monomers, then dimer concentrations dominate at later stages (Momiji and Monk, 2008), which results in slower decay at later time points. (2) If the proteins triggering degradation of Her7 are also short-lived, cycloheximide treatment could prevent the replenishment of the degradation complex. These options will be investigated in future studies. The measured value has to be treated as an upper limit for the half-life of Her7 as: (1) we overexpressed the HA-Her7 protein above its normal levels, and (2) the short HA-tag might have mildly stabilized Her7 protein. Nevertheless, the measured half-life corroborates the prediction of our model.

Oscillators are one of the most common network modules that regulate the timing of biological processes. Mechanistic understanding of the segmentation clock will clearly shed light into the comprehension of oscillators in other biological systems. Our study provides a comprehensive tissue-level stochastic modeling of the system explaining the amplitude, period and synchronization dynamics of the zebrafish somite segmentation clock. Although direct proof for the her autoinhibitory feedback loop acting as the pacemaker of the zebrafish segmentation clock is still lacking, our multicellular stochastic model suggests that the transcriptional feedback-based model could be sufficient to describe the inner workings of the clock. As multiple oscillations are running in parallel in amniotes, it is not clear whether the segmentation clock is really based on a simple transcriptional feedback loop in either zebrafish or amniotes (reviewed by Ozbudak and Pourquié, 2008). Thus, direct experimental tests are necessary to assess the accuracy of transcriptional feedback-based models (such as Goldbeter and Pourquié, 2008) in all vertebrate model organisms.

We thank members of the Ay and Ozbudak laboratories for their comments. We are grateful to Julian Lewis and François Giudicelli for providing the Tg(hsp70l:HA-her7) transgenic line and Spartak Kalinin for fish care. We thank Anja Hanisch, Dan Schult, Ken Belanger, Stuart Newman and Julian Lewis for critical reading of the manuscript.