## Abstract

The transcriptional repressor Hes1, a basic helix-loop-helix family protein, periodically changes its expression in the presomitic mesoderm. Its periodic pattern of expression is retained in a number of cultured murine cell lines. In this paper, we introduce an extended mathematical model for Hes1 oscillatory expression that includes regulation of Hes1 transcription by Drosophila Groucho (Gro) or its vertebrate counterpart, the transducine-like enhancer of split/Groucho-related gene product 1 (TLE1). Gro/TLE1 is a necessary corepressor required by a number of DNA-binding transcriptional repressors, including Hes1. Models of direct repression via Hes1 typically display an expression overshoot after transcription initiation which is not seen in the experimental data. However, numerical simulation and theoretical predictions of our model show that the cofactor Gro/TLE1 reduces the overshoot and is thus necessary for a rapid and finely tuned response of Hes1 to activation signals. Further, from detailed linear stability and numerical bifurcation analysis and simulations, we conclude that the cooperativity coefficient (*h*) for Hes1 self-repression should be large (i.e. *h*≥4). Finally, we introduce the characteristic turnaround duration, and show that for our model the duration of the repression loop is between 40 and 60 min.

## 1. Introduction

Recent experiments on cultured mammalian cell lines have shown oscillatory dynamics in the expression of at least three transcriptional factors: Hes1, p53 and NF-κB (Lev Bar-Or *et al*. 2000; Hirata *et al*. 2002; Hoffmann *et al*. 2002). In this paper, we deal with one of them, Hes1. We propose an extended model for the regulation of Hes1 expression that explains its oscillatory nature, and look for parameter requirements in different versions of the model.

Somites are transient embryonic structures that are formed through segmentation of the presomitic mesoderm (PSM) in a highly regulated process called somitogenesis. They are the origin of the skeletal muscles of the body as well as the axial skeleton and the dermis of the back (Christ & Ordahl 1995). The partitioning of the vertebrate body into a repetitive series of somites requires the spatially and temporally coordinated behaviour of the cells in PSM. This is achieved by an oscillatory mechanism generating periodic waves of gene expression (order of 2 h) known as segmentation clock (Schnell *et al*. 2002).

The genes that are periodically expressed belong to the HES family (for Hairy Enhancer of Split; Gao *et al*. 2001). They act by negatively regulating transcription of tissue-specific transcription factors (Ohsako *et al*. 1994). The transcription factor Hes1 is essential to neurogenesis, myogenesis, haematopoiesis and sex determination (Proush *et al*. 1994; Ishibashi *et al*. 1995). Hirata *et al*. (2002) showed that a single serum treatment induces an oscillation of Hes1 mRNA and protein concentrations with a 2 h period in a variety of cultured cells, such as myoblasts, fibroblasts, neuroblastoma cells and teratocarcinoma. They also showed that Hes1 acts as its own repressor through a negative feedback loop. To account for the appearance of these oscillations, they proposed the involvement of a third factor in addition to Hes1 protein and Hes1 mRNA, because three species are necessary to induce oscillations in a negative feedback system (Griffith 1968). However, the introduction of an explicit delay of 15–20 min in a basic two-species model with negative feedback is also sufficient to explain the oscillations (Jensen *et al*. 2003; Monk 2003). The delay accounts for the time needed for transcription, translation and formation of a complex in the nucleus to start repression.

Drosophila Groucho (Gro) and its vertebrate counterparts, the transducine-like enhancer of split/Groucho-related gene products 1–4 (TLEs1–4), lack DNA binding ability but can functionally associate with a number of DNA-binding proteins, including Hes1. Interaction with Gro/TLE1 at the nuclear matrix is necessary for transcriptional repression by Hes1 (McLarren *et al*. 2001). Hes1 itself activates hyperphosphorylation of Gro/TLEs bound to it. This correlates with association of Gro/TLE1 to a nuclear matrix, suggesting that chromatin remodelling might be a mechanism for this transcriptional repression (Nuthall *et al*. 2002). We present here a model that includes Gro/TLE1 in the regulation of Hes1-mediated repression of transcription.

Ordinary differential equations (ODEs) with negative feedback loops need to be at least three-dimensional to sustain oscillations. A requirement for these oscillations is a high cooperativity coefficient, i.e. a large feedback gain is necessary for the onset of oscillations (Tyson & Othmer 1978). One such system, the Goodwin model (Goodwin 1965), has been extensively used in circadian clock modelling. However, this model requires a cooperativity coefficient of at least 9 to induce oscillations, which is large and unrealistic unless a cascade mechanism is involved (Ferrell 1996). On the other hand, the introduction of a time delay in a feedback loop allows sustained oscillations even in a one-dimensional delay differential equation with a relatively low cooperativity coefficient.

Time-series of Hes1 mRNA and protein levels in various cell cultures were used as a starting point for modelling Hes1 repression loop. This paper considers two different models of self-repression kinetics, and uncovers conserved features among them. Section 2 presents an analysis of the two-dimensional system proposed by Jensen *et al*. (2003). In §3, we develop an extension of the Jensen *et al*. model considering Gro/TLE1 interaction with the Hes1 regulatory pathway. We then discuss the effects of Gro/TLE1 on the properties of oscillatory Hes1 expression after initiation by a serum shock. Section 4 considers general features of the two models that allow the estimation of parameter values. We find quantities that are relatively invariant with respect to different models. We conclude that although interaction of Hes1 with the corepressor Gro/TLE1 is not necessary for oscillatory expression of Hes1, it does allow a much more realistic activation curve that is consistent with the experimental data.

## 2. Two-dimensional Hes1–mRNA repression model

Standard mathematical analysis shows that two-component models with a negative feedback cannot have stable self-sustained oscillations. Two possibilities are then offered for a simple extension of the model to achieve oscillatory behaviour.

First, one can invoke the existence of a third dependent variable. This is the route taken by Hirata *et al*. (2002) who have considered a three-dimensional system, in which a third (unknown) species *X* actively degrades Hes1, while Hes1 also acts as a repressor for *X*. However, the analysis of this model is not further considered in this paper. Second, an explicit delay (Jensen *et al*. 2003; Monk 2003; Lewis 2003) is considered, which has its origin in the underlying biology. In either case the Hill coefficient *h*, which accounts for the cooperativity of repression activity of agents involved in the inhibition of Hes1 transcription, is crucial. Although these studies dealt with parameter estimation and stability analysis, the cooperativity coefficient *h* was taken as a fixed parameter with a value between 2 and 4. The complexity of transcription regulation does not allow for direct evaluation of the cooperativity coefficient, so we analysed a self-repression to find plausible ranges of the value of *h*.

A two-compartment model for Hes1–mRNA self-repression with delay can be written as (Jensen *et al*. 2003; Monk 2003; Lewis 2003)(2.1)(2.2)Equations (2.1) and (2.2) will be referred to as Model A. The first equation describes the cellular concentration of Hes1 mRNA () at time *t* and the second equation describes the cellular concentration of Hes1 protein () at time *t*. The respective degradation rates (*δ* of and *α* of ) and translation rate (*β*) all have units of min^{−1}. The nonlinearity comes from the transcriptional repression activity of Hes1. The parameters in the feedback loop are the maximum mRNA transcription rate (conc min^{−1}), a DNA dissociation constant *k* (same units as ) and a Hill coefficient *h* representing the degree of cooperativity between different factors intervening in the repression of *Hes1* gene by Hes1 protein. The delay *τ* accounts for the delay due to Hes1 modification, complex formation and nucleocytoplasmic molecular transport.

Table 1 gives the values of parameters used for the analysis and simulations. The experimentally measured parameters are the decay rates *α* and *δ*, and the period of oscillation . The deactivation rate *σ* was assumed to be of the same order as the degradation rates. We have performed a numerical sensitivity analysis by varying every parameter from 0.1 to 10 times their default value. We looked at the presence of oscillations and the period of oscillation. We found that the oscillations were present for most of the parameter values and that the period changed by less than 15% with respect to , , *k*, and *β* (results not shown). Degradation rates had a marked influence on the behaviour of the system, but they have been measured accurately (Hirata *et al*. 2002). Therefore, the following analysis mainly focuses on the effect of the Hill coefficients *h* and *u* and the explicit delay *τ*. The details of the computations and analysis are given in appendices A–C.

The bifurcation diagram with respect to *h* (figure 1*a*) indicates the two possible behaviours of Hes1 activity: damping to a steady state or sustained oscillations, with an increasing period (figure 1*c*). Experimental levels of Hes1 mRNA and protein oscillate with a relatively constant amplitude and a 2 h period (figure 1*b*).

In this model, a stability analysis (cf. appendix A) leads to analytical formulae for local stability connecting the Hill coefficient, Hopf (oscillation) period and rate constants. Namely, the steady state is unstable and a limit cycle is obtained whenever *h* is larger than(2.3)Note that this critical value of *h* does not depend on the delay *τ*. Rather, we use the Hopf period as a parameter and let *τ* depend on . This allows us to place bounds on *h* and *τ*. Assuming a period , we obtain a lower bound on the cooperativity coefficient . For a given cooperativity coefficient *h* and Hopf period , the associated critical delay can be expressed as(2.4)where *ω* is the oscillatory angular frequency, . The critical delay defines the boundary between stable () and oscillatory () mRNA and protein levels. This leads to a value of the maximal , .

When Model A has no explicit delay (*τ*=0), its characteristic time-scale (CTS), defined as with *λ* the eigenvalue of Model A, is determined by the decay rates *α* and *δ* (Murray 1993; Strogatz 1994). Explicitly, . Added to , it gives a characteristic turnaround duration (CTD) of . The turnaround duration is the time required to have an effective repression after transcription initiation (see appendix B for a derivation of these results).

The value constitutes an upper bound on the delay and is a lower bound on the cooperativity coefficient. To have large amplitude oscillations (i.e. to move away from the bifurcation point), it is necessary to have a somewhat smaller *τ* and/or larger *h*.

## 3. Gro/TLE1-mediated repression allows tuned response

### (a) Gro/TLE1–Hes1 repression model

We now consider the influence of an additional factor known to be involved in the Hes1 repression loop, namely Gro/TLE1. Protein Gro/TLE1 is activated through Hes1-induced hyper-phosphorylation. This activation is described by a Hill function that is a monotonically increasing function of Hes1 with Hill coefficient *u*. Moreover, the active form of Gro/TLE1 forms a complex with Hes1, denoted *GroH*, to mediate repression through a negative feedback loop (cf. figure 2). The variable *GroH* represents the repression complex of hyper-phosphorylated Gro/TLE1 with Hes1 protein. The associated equations are written as(3.1)(3.2)These equations constitute what we call Model B. As with the other model, local stability considerations establish relationships between the Hill coefficient, period, kinetic constants and delay. Following the same analytical development as for Model A, we found that the quantity playing the role of the Hill coefficient was now the product of *h* and *u*. Denoting this quantity total cooperativity, , we obtain an expression similar to equation (2.3):(3.3)Again, we can calculate the associated critical delay ,(3.4)As in Model A, a threshold value of the total cooperativity also defines the boundary between stable steady state () and limit cycle (). The introduction of another nonlinearity does not change the linear stability analysis performed for Model A. It should be noted, however, that this analysis is an approximation and a more detailed numerical analysis shows that the Hill coefficient threshold is slightly higher than predicted (2.5 instead of 2.0 for fixed ). The CTD has, as for Model A, a value of CTD=53 min.

Figure 3 shows the stability boundary in space for Models A and B (figure 3*a*), and the corresponding Hopf period (figure 3*b*). The region of oscillatory expression is above and to the right of the curves. For Model B, only the Hill coefficient *h* is shown. With *u*=2, the total cooperativity coefficient is doubled, leading to almost the same result as for Model B. Thus, for a 120 min period oscillation, the corresponding values of *h* and *τ* can be uniquely determined. For Model B, the Hill coefficient required for 120 min period oscillations is greater than 2.5 and the delay smaller than . If an interval between 90 and 150 min for the period is considered, a range for *h* and *τ* can be defined as:

for Model A, and ;

for Model B, , and .

These results show that introducing a new mechanism does not affect the basic properties of these models with respect to the Hopf bifurcation. Interestingly, the model considered by Hirata *et al*. (2002) consists of three ODEs with two feedback loops, each with a Hill coefficient of 2 that would give a total cooperativity coefficient of 4. Our comparison of these models reveals general properties of the Hes1 oscillation that are not yet measured experimentally. In particular, strong cooperativity of Hes1 repression (*h* between 3 and 8) and a turnaround duration between 40 and 60 min are predicted.

### (b) The second nonlinearity increases adaptativity

The stability results of the previous sections show that oscillation onset of the systems depends on model independent quantities, such as degradation rates, cooperativity coefficients (or the product thereof) and the CTD, including delays due to auxiliary variables. The CTD takes into account both the time-scale introduced by the dynamics of the non-delayed system and the explicit delay *τ*. However, the behaviour of the solution away from the steady state can be quite different from model to model.

Experimental data on Hes1 expression level in cultured cells after a serum shock show no sign of overshoot in the first cycle of mRNA transcription and Hes1 synthesis (Hirata *et al*. 2002). On the contrary, the expression levels of both mRNA and Hes1 rapidly settle down to an oscillatory regime of approximately constant amplitude. Systems with only one nonlinearity often display a large overshoot before solutions converge to their attractor. A second nonlinearity is then needed to fine-tune the dynamics of the systems when away from equilibrium. In appendix C, using Model A, we compute a lower bound for the value of the Hes1 expression overshoot with respect to the steady-state value.

From the results of a full numerical simulation of Models A and B, we compared the initiation of Hes1 synthesis after a serum shock (figure 4). The serum shock triggers Hes1 transcription and is modelled by setting , and initial levels to values close to zero at time 0. In Model A, there is an overshoot due to the lack of repression mechanisms in the first minutes. In the case of Gro/TLE1-mediated repression (Model B), it is assumed that Gro/TLE1 is already active at a low level, and the fast activation rate of Gro/TLE1 supports the rapid convergence of Hes1 expression to the stable oscillatory pattern.

It is possible to give a formal bound for the overshoot. The following calculation shows that the large overshoot (figure 4) is a characteristic feature of Model A. We assume that before initiation the *Hes1* gene is practically silent, so the initial conditions for and are set to zero. Moreover, to simplify the computations, we set *α*=*δ*. To obtain a simple expression for the overshoot, we set the transcriptional delay *τ* equal to the half-lives of Hes1 mRNA and protein, so . Calculations of a bound of the overshoot in Hes1 lead to (see appendix C). With a steady state of , the ratio between the maximal protein level and the steady-state level is(3.5)The ratio of maximal transcription to degradation rates, , is large. For instance, gives an overshoot about five times the steady-state value. It is likely that the ratio is even higher since estimates of 11 initiations per second have been reported (Bolouri & Davidson 2003). It is possible to reduce the overshoot considerably by setting some parameters to unrealistic values. For example, a translation rate *β* of one translation every 33 min (0.03 min^{−1}) leads to a smaller overshoot, comparable to Model B simulations. It is likely, however, that *β* is much larger, ranging from 1.0 to 100 min^{−1} (Ghaemmaghami *et al*. 2003; Bolouri & Davidson 2003). Even though *β* has no effect on the stability properties of these models, initially it is of crucial importance after serum shock.

In principle, one could replace *β* in Model A by a nonlinear function without having to take Gro/TLE1–Hes1 interactions into consideration. However, the introduction of this nonlinearity is not motivated by any experimental data and would be artificial.

Compared to Model A, Model B shows a more realistic activation curve. This is due to two factors. First, since Gro/TLE1 is a general corepressor, there may exist a baseline level of phosphorylated protein in the cell before initiation of Hes1 synthesis. This would have an attenuating effect on Hes1 transcription. Second, numerical simulations show that in Model B a smaller delay *τ* is needed to generate 120 min period oscillations due to the introduction of a second nonlinearity. Consequently, the speed of the response after activation is increased.

## 4. Discussion

Numerical simulations and analytical results from two models of Hes1 self-regulation show interesting model-independent features. Encouragingly, the most important characteristics of the models all lie in the same range. More precisely, to be consistent with the experimental data the cooperativity coefficient associated with Hes1 repression must be higher than 2, and is likely to be around 4. Further, the explicit delay included in the model should range from 10 to 30 min. These results are in agreement with previous studies (Hirata *et al*. 2002; Jensen *et al*. 2003; Monk 2003; Lewis 2003). Given the complex nature of transcription regulation, it is likely that this parameter plays an important role in the control of protein expression.

A central quantity introduced in this study is the CTD, defined as the CTS+the delay, as a measure of the lag between transcription initiation and repression. In both model versions, a CTD of about 40–60 min is predicted. The model-independent estimation of is a central result based only on partial knowledge of kinetic data. A sufficiently large delay is needed to induce oscillations, so the CTS needs to be small. Since the CTS is , only protein and mRNA with short half-lives can display these short oscillation periods.

Using linear stability analysis, we estimated the values of the cooperativity coefficients and delays from measured kinetic data, such as Hes1 protein and mRNA degradation rates. Constrained by the period of oscillation, it is possible to define a lower bound for the cooperativity coefficient and an upper bound for the delay. For example, a cooperativity coefficient lower than 4 or a delay longer than 20 min always yield periods longer than 120 min. If the system is near the Hopf bifurcation, these values should be good estimates.

Interaction of Hes1 with Gro/TLE1 is necessary for transcriptional repression by Hes1 (McLarren *et al*. 2001). The fact that Hes1 itself activates hyperphosphorylation of Gro/TLEs suggests a dynamic link between the two proteins. Introduction of Gro/TLE-mediated repression in a model allows a faster physiological adaptation after a serum shock and/or Hes1 induction from the Notch pathway. In feedback loop systems with a single nonlinearity, a large overshoot is seen after initiation, resulting in an initial strong response followed by small amplitude oscillations. When the corepressor Gro/TLE1 is introduced, the overshoot is greatly reduced, and steady oscillatory protein expression levels are observed. Gro/TLE1 is known to be expressed in a variety of cell lines and acts as a general corepressor (Jimenez *et al*. 1997). We suggest that the kinetic role of Gro/TLE1 is to fine-tune expression levels after initiation of protein synthesis.

## Editors' note

Please see also related communications in this focussed issue by Benson *et al*. (2006) and Terashima *et al*. (2006).

## Acknowledgments

This work was supported by MITACS (Canada) and the Natural Sciences and Engineering Research Council (NSERC grant OGP-0036920, Canada). We especially thank Drs M. Santillán, S. Stifani and R. Kageyama for helpful discussions and comments. This work was initiated while B.Č. was visiting the Center for Nonlinear Dynamics, McGill University. We acknowledge support from the German grant agencies BMBF and DFG.

## Footnotes

One contribution of 13 to a Theme Issue ‘Biomathematical modelling I’.

- © 2006 The Royal Society