Modelling the checkpoint response to telomere uncapping in budding yeast

C.J Proctor, D.A Lydall, R.J Boys, C.S Gillespie, D.P Shanley, D.J Wilkinson, T.B.L Kirkwood

Abstract

One of the DNA damage-response mechanisms in budding yeast is temporary cell-cycle arrest while DNA repair takes place. The DNA damage response requires the coordinated interaction between DNA repair and checkpoint pathways. Telomeres of budding yeast are capped by the Cdc13 complex. In the temperature-sensitive cdc13-1 strain, telomeres are unprotected over a specific temperature range leading to activation of the DNA damage response and subsequently cell-cycle arrest. Inactivation of cdc13-1 results in the generation of long regions of single-stranded DNA (ssDNA) and is affected by the activity of various checkpoint proteins and nucleases.

This paper describes a mathematical model of how uncapped telomeres in budding yeast initiate the checkpoint pathway leading to cell-cycle arrest. The model was encoded in the Systems Biology Markup Language (SBML) and simulated using the stochastic simulation system Biology of Ageing e-Science Integration and Simulation (BASIS). Each simulation follows the time course of one mother cell keeping track of the number of cell divisions, the level of activity of each of the checkpoint proteins, the activity of nucleases and the amount of ssDNA generated. The model can be used to carry out a variety of in silico experiments in which different genes are knocked out and the results of simulation are compared to experimental data. Possible extensions to the model are also discussed.

Keywords:

1. Introduction

Telomeres are repetitive sequences of DNA situated at the ends of linear chromosomes. They require protection from being recognized as double-strand breaks to prevent activation of the DNA damage response. In budding yeast, the telomere-binding protein Cdc13 provides a protective cap at the ends. In the absence of Cdc13, checkpoint proteins bind to the telomeres, resulting in cell-cycle arrest (Weinert & Hartwell 1988; Lydall 2003), accumulation of single-stranded DNA (ssDNA) on the 3′ strands (Garvik et al. 1995) and recombination (Grandin et al. 2001a). A mutant form of Cdc13, cdc13-1, has a temperature-dependant response and is only able to cap telomeres at a range of temperatures, with the maximum ‘permissive temperature’ being 25°C (Garvik et al. 1995). At higher temperatures, cdc13-1 is unable to function, the checkpoint pathway is activated and the cells stop growing. Hence, a temperature in this range is referred to as a ‘restrictive temperature’. A useful way to study the checkpoint response in budding yeast is to introduce the cdc13-1 mutation and to initially culture the cells at a permissive temperature (e.g. 23°C) and then switch to a restrictive temperature (e.g. 36°C). Throughout this paper, we will use the terms permissive temperature and restrictive temperature to refer to cdc13-1 strains being cultured at 23 and 36°C, respectively. More detail of the experimental framework can be found in Zubko et al. (2004).

The end of a telomere in budding yeast consists of a very short region of ssDNA (approx. 12 nt) which is termed an ‘overhang’ (Larrivee et al. 2004). This overhang is usually protected by telomere-binding proteins to prevent it from being seen as DNA damage, which would otherwise activate a DNA damage response. At the restrictive temperature, cdc13-1 cells possess a small amount of DNA damage, as a result of overhangs becoming exposed, which activates the Rad9- and Rad24-dependent checkpoint pathway. This initial damage is amplified to larger single-stranded regions near telomeres. The generation of ssDNA is a necessary requirement for the full activation of the checkpoint response. There is evidence that the Mrx complex, Exo1 and the Rad24 group regulate or encode nucleases, which are responsible for generating ssDNA at telomeres (Lydall 2003). Exo1 regulates ssDNA levels when telomere capping is defective. For example, Exo1 is essential for ssDNA generation in yku70Δ mutants and contributes to ssDNA in cdc13-1 mutants (Maringele & Lydall 2002). In cdc13-1 rad9Δ mutants, the ssDNA accumulates more rapidly than in cdc13-1 mutants, whereas in cdc13-1 rad24Δ mutants, the ssDNA accumulates much more slowly than in cdc13-1 cells (Lydall & Weinert 1995). It has been found that ssDNA levels are high in cdc13-1 rad9Δ exo1Δ strains, low in cdc13-1 exo1Δ strains and low in cdc13-1 rad53Δ exo1Δ and cdc13-1 mec1Δ exo1Δ strains (Jia et al. 2004; Zubko et al. 2004). This suggests that Rad9 inhibits an unidentified nuclease called ExoX, and that this inhibition does not require Mec1 or Rad53 (see figure 1).

Figure 1

A model for the interaction between checkpoint pathways and nucleases at cdc13-1-induced damage (Reproduced with permission from Jia et al. (2004), Copyright Genetics Society of America.). Rad9 inhibits ssDNA production at telomeres of cdc13-1 strains by two routes, a and b. One route, a, depends on Mec1, Rad9 and Rad53 and targets Exo1. The other route, b, is dependent on Rad9, independent of Mec1 and Rad53 and targets ExoX.

In order to understand precisely how this important yet complex system of pathways is regulated, there is a clear role for mathematical models. These can help both to reveal gaps and inconsistencies in current knowledge and also to conduct in silico experiments (simulations) that can help in the efficient planning of wet experimental studies. A model for the interaction of checkpoint pathways and nucleases at cdc13-1-induced damage has been proposed by Jia et al. (2004). The model from their paper has been reproduced in figure 1. The Rad17, Mec3 and Dcd1 complex is a hetero-trimeric proliforating cell nuclear antigen (PCNA)-type ring, which is loaded onto damaged telomeres by another complex consisting of Rad24 and Rfc2–5 subunits, which together form a clamp–clamp loader system (Majka & Burgers 2003). Once loaded, the Rad17 complex activates ExoX, which generates ssDNA at the telomeres. Mec1 and Ddc2 bind ssDNA independent of the Rad17 complex. Mec1 is essential for phosphorylation and activation of Rad9, which in turn activates two parallel pathways of cell-cycle arrest. One pathway is via the kinase Chk1 and the other pathway is via the kinase Rad53 (Gardner et al. 1999; Sanchez et al. 1999). Rad9 also inhibits the activity of ExoX (Zubko et al. 2004). Exo1 activity is not dependent on Rad24 or Rad17, but its activity is inhibited by Rad53. It has been proposed that Rad9 mediates interactions between the upstream kinase, Mec1 and two parallel downstream kinases, Rad53 and Chk1 (Gilbert et al. 2001; Schwartz et al. 2002; Blankley & Lydall 2004).

In a study of yeast checkpoint genes, it was found that there was a correlation between the amount of ssDNA and cell death in several single and double mutant strains (Lydall & Weinert 1995). However, more recently it has been found that mec1Δ, rad53Δ and exo1Δ mutations each suppress the rapid loss in viability of cdc13-1 rad9Δ mutants, but ssDNA still accumulates (Jia et al. 2004). They suggest that ssDNA is only cytotoxic if Mec1, Rad53 and Exo1 convert it into a lethal lesion.

The biochemical mechanisms underlying the DNA damage response are complex and many hypotheses have been put forward to suggest how DNA repair and checkpoint pathways interact. This has motivated us to develop a mathematical model of the system. To build the model, it was necessary to specify precisely each element of the network and its interaction with other elements (see §2.1 for details).

2. The model

2.1 Elements of the model

The paper describes a core model representing the wild-type from which a number of variants were created to represent a number of experimental knockouts. The set of models, all contain the same network, but the variants have different sets of parameter values or initial conditions. The differences are typically small, but lead to significant variance in phenotype. We use the term ‘species’ to denote a biochemical entity in the model and interactions between species are referred to as reactions. These terms will be used throughout this paper. The first step in building the model is to define all the species and specify their initial amounts. A complete list for the wild-type variant is given in table 1 and each of the species and all reactions in the model are described later in detail. A complete list of the reactions and their parameters (for the wild-type variant) is given in table 2. We use mass action stochastic kinetics for the rate laws (see Wilkinson 2006 for further details). We also use event structures in our model, which allows species amounts and parameter values to be changed once a particular condition is true. These are detailed in the text and summarized in table 3. Figure 2 shows the reaction network in the model.

Figure 2

Network diagram showing the checkpoint response to uncapped telomeres. (a) Activation of ExoX and Exo1. ExoX requires Rad24 and Rad17 binding for its activation. Exo1 is activated independent of Rad24 and Rad17, although it may also act on telomeres bound by Rad17 and be activated in a Rad24-dependent manner. (b) Binding of single-stranded DNA (ssDNA) by RPA. Each molecule of RPA requires three units of ssDNA to bind. (c) Activation of checkpoint response via Rad53/Dun1 pathway. Activation of this pathway leads to inhibition of nuclease activity. Dashed line indicates an event, where a threshold level of Mec1RPAssDNA activates a kinase (Rad9Kin) which activates Rad9. (d) Activation of checkpoint response via Chk1 pathway does not affect Exo1. (e) Recovery can take place during S phase or G2/M arrest after single-stranded DNA has been removed. The dashed line indicates an event. When the level of ssDNA is equal to 0, the dummy species ‘recovery’ is set to 1 and the recapping reaction can then occur. See §2.1.6 for more details.

View this table:
Table 1

Initial conditions for wild-type model.

View this table:
Table 2

Reactions for wild-type model.

View this table:
Table 3

Summary of events used in all models.

2.1.1 Capped telomeres, uncapped telomeres and Cdc13

Budding yeast have 16 chromosomes and hence 32 telomeres in G1, but during S phase, the DNA is replicated resulting in 64 telomeres. To avoid unnecessary complexity, we have omitted detail of DNA replication from this model and assumed that the number of telomeres is constant. The main purpose of this model is to examine the pathways that lead to G2/M arrest and so we have chosen 64 telomeres. Telomeres can be in either a capped or uncapped state and hence we have two species to represent telomeres, which we call Utelo and Ctelo for uncapped and capped telomeres, respectively. We assume that the ssDNA overhangs at the telomere ends are capped by Cdc13. The dissociation rate (the time required for half the protein–DNA complex to dissociate) is approximately 30 min (Lin et al. 2001). Cdc13 has a strong binding affinity to single-stranded telomeric DNA with a dissociation constant of ca 10−7 M (Lin et al. 2001). Hence, there is equilibrium between uncapped telomeres (Utelo) and those bound by Cdc13.Embedded Image

To model this equilibrium, we set up two reactions,Embedded Imagewhere Ctelo represents the pool of telomeres bound by Cdc13.

The first reaction represents capping and the second reaction uncapping. The capping reaction is a second-order reaction since there are two reactants and is given by k1[#Utelo][#Cdc13], where ‘#’ represents the number of molecules. The uncapping reaction is a first-order reaction and is given by k2[#Ctelo]. The value of k2 can be calculated using information on the time taken for half of the bound Cdc13 to dissociate from the telomere given earlier. Since the biological half-life is given by t1/2=−ln(0.5)/degradation rate, it follows that Embedded Image. The dissociation constant kd is 10−7 M=110 molecules per nucleus, assuming that the volume of a yeast nucleus is ca 1.8 μm3 (1.8×10−9 μl; http://biochemie.web.med.uni-muenchen.de/Yeast_Biol/). Since kd =k2/k1, then Embedded Image. We assume that all the telomeres are in a capped state initially, but there is a dynamic process of uncapping and recapping, hence at any point in time, there is a very low number of uncapped telomeres. Telomeres remain in an uncapped state only very transiently, and signalling proteins do not bind while they are in this transient state.

2.1.2 Binding of uncapped telomeres: Rad17 and Rad24

If a telomere remains uncapped, then the Rad17 complex is loaded onto the telomere end by the Rad24 clamp loader complex. This loading requires ATP and can be represented by the biochemical reactionEmbedded Imagewhere Rad17Utelo represents the binding of Rad17 to the telomere end. For simplicity, we have omitted details of Mec3 and Ddc1 in the Rad17 complex and omitted the Rfc2–5 subunits in the Rad24 complex, since they are unimportant for this model. Since there are four reactants in the Rad17-binding reaction, the stochastic rate law is a fourth-order reaction. Care needs to be taken when considering higher order reactions. For example, it is not realistic that the reaction rate would increase indefinitely with an increase in ATP levels; doubling the level of ATP in the cell would not necessarily double the reaction rate. Instead, we would expect that the reaction rate would increase asymptotically with ATP levels. Therefore, we use the rate law k3[#Utelo][#Rad17][#Rad24][#ATP]/(5000+[#ATP]).

The term in the denominator ensures that the reaction rate increases approximately linearly with ATP when ATP levels are low (less than ca 2000), but for high levels of ATP (greater than 5000), the increase in the reaction rate starts to level off with increasing ATP levels. The value for k3 was chosen by making an initial estimate and then simulating the model and comparing output to the experimental data. Further details of how we chose parameter values can be found in §2.3.

When Cdc13 is present, this reaction is very unlikely to occur as the binding strength of Rad17 compared to Cdc13 is extremely low. However, in the absence of Cdc13, this reaction will eventually take place and an irreversible set of reactions is initiated leading to the generation of ssDNA and cell-cycle arrest unless the pathway is interrupted downstream. We chose a very low value for k3, to ensure that the probability of a telomere, which is only transiently uncapped, binding to Rad17 is very low.

2.1.3 Nuclease activity: ExoX, Exo1 and ssDNA

ExoX is recruited to telomeres bound by Rad17, and degrades one of the strands to leave a long stretch of ssDNA. This reaction can only occur if ExoX is in its active state. The binding of telomeres by Rad17 (which in turn requires Rad24) is required for ExoX activation. This can be represented simply by the following reactions:Embedded Imagewhere ExoXI represents the inactive form of ExoX and ExoXA represents activated ExoX. We assume that each time this reaction takes place, 10 nt ssDNA are produced and we set one unit of ssDNA equal to 10 nt. The reaction can continue to take place until ExoX is inactivated. It has been observed that ca 8 kb ssDNA are generated in 1 h in cdc13-1 strains (Jia et al. 2004; Zubko et al. 2004). Exo1 degrades not only unprotected telomeres, independent of ExoX and Rad17, but also telomeres bound by Rad17. Exo1 activation does not require Rad24. Since cdc13-1 rad24Δ mutants have low levels of ssDNA (Lydall & Weinert 1995), we assume that activation of Exo1 is much lower than that of ExoX. There are very low levels of ssDNA in the cdc13-1 rad24Δ strain, which must be due to Rad24-independent activation and subsequent activity of Exo1. However, the difference between the amounts of ssDNA in cdc13-1 exo1Δ compared to cdc13-1 mutants is large (Zubko et al. 2004), suggesting that Exo1 activity is higher when Rad24 is present. Therefore, we also include an additional reaction of Rad24-dependent activation of Exo1, although this was not in the model proposed by Jia et al. (2004) (see figure 1). The reactions for nuclease activity are as follows:Embedded Image

2.1.4 Binding of ssDNA: replication protein A and Mec1

Replication protein A (RPA) binds to ssDNA with every molecule of RPA binding to 30 nt ssDNA (Zou & Elledge 2003). As each unit of ssDNA is equivalent to 10 nt, we assume that one molecule of RPA binds to three units of ssDNA. Mec1 in complex with Ddc2 then binds to RPA (Zou & Elledge 2003). Thus, binding of ssDNA can be represented by the following two reactions:Embedded Image

Since three units of ssDNA are required for the RPA-binding reaction, the mass action stochastic kinetic rate law is given by k8a[#RPA][#ssDNA][#ssDNA-1][#ssDNA-2]/6. However, this rate law is only valid if the ssDNA units are independent. This is certainly not the case, and the rate law should be a linear and not a cubic function of the total amount of ssDNA. Therefore, we model RPA binding as a series of steps. To do this, we introduce two dummy species RPAssDNA1 and RPAssDNA2 for the intermediate steps. In the first step, RPA binds to one unit of ssDNA to form RPAssDNA1. The rate of this reaction depends linearly on the total amount of ssDNA. In the next step, RPAssDNA1 binds to the second unit of ssDNA to form RPAssDNA2 and finally RPAssDNA2 binds to the third unit of ssDNA to form RPAssDNA. Once the first step has happened, the next two steps follow quickly by setting the rate constants to a high value (k8b=k8c=100). The full set of reactions for RPA and Mec1 binding areEmbedded Image

2.1.5 Checkpoint activation: Rad9, Rad53, Chk1 and Dun1

A very large threshold of ca 10 kb ssDNA is required for keeping checkpoints active in order to maintain cell-cycle arrest (Vaze et al. 2002; Zubko et al. 2004; Maringele & Lydall 2005). To prevent activation of Rad9 when only low levels of ssDNA are present, we modelled this activation as a two-step process. First, we used an event structure in the model to activate a kinase that activates Rad9 (called ‘Rad9Kin’), if the total amount of ssDNA within the cell exceeded 8 kb. The rate for the Rad9-activation reaction is equal to k9 [#Rad9I][# Rad9Kin] and will be equal to 0 as long as Rad9Kin=0. Once Rad9Kin is activated (so that Rad9Kin=1), the Rad9-activation reaction can proceed. The event and reaction representing Rad9 activation are as follows:Embedded Image

Activated Rad9 activates both Rad53 and Chk1. Rad53 in turn, inhibits the activity of Exo1 and activates Dun1. Rad9 also inhibits ExoX. This inhibition does not require Rad9 to be activated by Mec1 and so it is represented by two reactions, one with Rad9I and the other with Rad9AEmbedded ImageActivation of Chk1 and Dun1 results in G2/M arrestEmbedded Imagewhere G2Mon and G2Moff are dummy species, which control whether the cell cycle can progress from G2 to M. If G2Mon=1 and G2Moff=0, then the above reactions can occur. If G2Mon=0 and G2Moff=1, then the reactions cannot take place (see appendix A for details of modelling the cell cycle).

2.1.6 Recovery

The ssDNA may be removed during S phase by DNA replication. It is also possible that ssDNA is removed during G2/M arrest. When all the ssDNAs have been removed, the uncapped telomeres can be recapped by Cdc13 and recovery takes place. After recovery, all the checkpoint proteins are turned off and cell division can take place. This procedure is modelled by a series of reactions and eventsEmbedded ImageIf Mec1RPAssDNA+RPAssDNA+ssDNA==0, then recovery=1 (where ‘recovery’ is a dummy species in the model).

If recovery=1, then the following reaction can proceed:Embedded ImageIf Rad17Utelo=0 (i.e. all telomeres have been recapped), then G2Moff is set to 0, G2Mon is set to 1, recovery is set to 0, and all other checkpoint proteins are reset to their initial conditions. A network diagram of the full system is shown in figure 2.

2.2 The cell cycle

We also include a series of reactions to represent the different stages of the cell cycle (see appendix A). We assume that ssDNA can be removed by DNA replication when the cell is in S phase and also to a lesser extent when cells are in G2/M arrest. Although it is unclear whether DNA replication can take place in G2/M arrest, it has been observed that the level of ssDNA does decrease in cells, which are arrested (Jia et al. 2004; Zubko et al. 2004). When budding yeast divides, the mother cell produces a small bud which then grows. When mitosis and cytokinesis are complete, the bud detaches to form the new daughter cell. This budding leaves a scar on the surface of the cell. Therefore, it is possible to count the number of cell divisions that a particular yeast cell has undergone by counting the number of budscars. We imitate this process in our model and keep track of the number of divisions by the introduction of a dummy species, budscar, which increases by one, every time the cell passes from M to G1 phase. The reaction and equation for the budscar species are shown in appendix A.

2.3 Initial conditions and parameter values

Before a model can be simulated, the initial amount of each species must be specified.

We obtained values for the majority of the initial amounts of proteins by consulting the yeast green fluorescent protein (GFP) fusion localization database (http://yeastgfp.ucsf.edu; Huh et al. 2003). We assume that initially all the checkpoint proteins are in their inactive state. In the G2 phase of the cell cycle, budding yeast contains 32 chromosomes and so have 64 telomeres, which we assume are all initially bound by Cdc13. We do not include detail of DNA replication in our model and so assume that the total number of telomeres is constant. Table 1 lists all the species, their respective initial amounts, the systematic names of all genes and the database terms so that each entity can be identified in public accessible bioinformatics databases.

The rate constants for each of the reactions are listed in table 2. Note that the value of k1 was set at 5×10−4, although we had calculated it to be 3.6×10−6 from published data (Lin et al. 2001), a 100-fold difference. This is because with our original value of k1, we found that the binding of Cdc13 to telomeres was not strong enough to keep them in a capped state (figure 3). Our model suggests that the binding of Cdc13 to telomeres must be stronger in vivo and it is probably the case that the binding affinity of Cdc13 is increased when it is in complex with Stn1 and Ten1 (Grandin et al. 2001b). At present, there are no experimental data to support this, but it seems a reasonable assumption. Therefore, we increased the value of k1 until we obtained results in which telomeres only become uncapped transiently as in figure 3 and chose k1=0.0005 as our default value. Figure 4 shows the model predictions for the kinetics of uncapping for the cdc13-1 strain when placed at the restrictive temperature using the default parameters. The only other kinetic data currently available for parameterizing the model are rates of single strand production in cdc13-1 mutants at high temperatures, where a rate of 8 kb h−1 has been observed (Jia et al. 2004; Zubko et al. 2004). The actual rates for the different nucleases are unknown, but they were chosen to reflect the suggested relative difference; therefore, we assumed that ExoX activity is 10 times faster than Exo1 activity.

Figure 3

Model predictions of the number of capped telomeres in a wild-type cell when k1=3.6×10−6 or k1=0.0005. The output is for one simulation over a period of 12 h.

Figure 4

Model predictions for the kinetics of uncapping in the cdc13-1 strain at the restrictive temperature (with the default parameters). The output is for one simulation over a period of 12 h.

We created an initial computer representation of the model using estimates for parameter values taken directly from evidence in the literature, where available. In some cases, parameter values started as initial best guesses. This initial model was then simulated and the results compared to experimental data. The parameter values were then adjusted to ensure that the model complied with that data. We are confident that the adjustments made to parameters were either to parameters for which we did not have specific data or that the data was derived from in vitro experiments, which do not necessarily correspond to the in vivo conditions of the culture experiments.

We coded the models in the systems biology markup language (SBML) level 2 (Hucka et al. 2003). SBML is a way of representing biochemical networks and is one of the standards used by the modelling community. It has been evolving since mid-2000 through the efforts of an international group of software developers and users. The models were developed using Mathsbml (Shapiro et al. 2004). The final model for the wild-type cell was encoded using SBML shorthand and then converted into full SBML (Wilkinson 2006). The model was then imported into the biology of ageing e-science integration and simulation (BASIS) system (www.basis.ncl.ac.uk) and the models for mutant strains were created by making simple modifications to the original model.

To simulate the model, we used a stochastic simulator based on the Gillespie algorithm (Gillespie 1977), since many of the processes being modelled occur at random times. This means that different outcomes may result from the same set of initial conditions. In addition, some of the species in our model must be considered as discrete units rather than continuous variables. For example, species representing the number of capped and uncapped telomeres must be a small integer and so it is not appropriate to use a deterministic simulator for this model. However, we also tried running simulations with a deterministic simulator to see if any parts of the model could be successfully modelled this way (see §3.3). The stochastic simulations were carried out on a Linux Beowulf cluster and the results were stored in a database. All the models and results in this paper are obtainable from the BASIS website (www.basis.ncl.ac.uk). For further details, see Gillespie et al. (2006a,b). The model for the wild-type strain, where cell death is included, is also available from the Biomodels database (ID:MODEL8679489165) at http://www.ebi.ac.uk/biomodels/ (Le Novère et al. 2006).

We did not carry out a full sensitivity analysis for the model parameters, since this would have been very time consuming due to the number of parameters and the fact that we used a stochastic simulator. However, many of the model parameters were chosen after many trials with different values and so it is possible to discuss the effects of changing these. The checkpoint response is sensitive to changes by an order of magnitude to the values of the parameters for the rate of capping (k1) and Rad17 binding to uncapped telomeres (k3). A 10-fold decrease in k1 or a 10-fold increase in k3 leads to an over-activation of the checkpoint response. The parameters involved in nuclease activity and removal of ssDNA were confined within fairly limited ranges in order to obtain rates of ssDNA production observed in cdc13-1 mutant strains. The parameters for cell-cycle arrest (k15 and k16) could be increased with less effect on the model predictions, but lowering the values led to an unsustained cell-cycle arrest, which could easily be overcome by Cdk activity. The rate of Rad9 activation (k9) was the most important parameter in the checkpoint pathway and affected not only how quickly the cell responded once the threshold level of ssDNA was reached, but also how quickly nuclease activity was inhibited. The other parameters in the checkpoint response had no significant effect on the model predictions for the rate of ssDNA production or the number of divisions before cell-cycle arrest or cell death. For example, a 10-fold increase or decrease on the rate of Mec1 binding (parameter k8d) did not change the rate of ssDNA production, the total level of ssDNA or the number of divisions before cell-cycle arrest.

3. Results

3.1 Wild-type cells

The model predicts that capped telomeres may become transiently uncapped, but are then recapped. At any time, there is usually only one uncapped telomere at the most, but on occasions, two or three telomeres may be uncapped simultaneously (figure 3). However, they do not remain uncapped long enough for Rad17 to bind, so very little ssDNA is produced (data not shown). Figure 5 shows the model predictions for the number of divisions obtained in three different simulation runs. Since we have used a stochastic simulator, each run produces a slightly different output. The cells divide ca 7–9 times in a period of 12 h, which gives an average division time of ca 1.5 h. This corresponds to the division time obtained when wild-type cells are grown at the restrictive temperature.

Figure 5

Model predictions for the number of divisions obtained in a wild-type cell. The output is for three simulations over a period of 12 h.

3.2 cdc13-1 mutant strains

3.2.1 Cell divisions

We used our model to simulate cdc13-1 mutant strains at the restrictive temperature by setting the initial amount of Cdc13=0, Ctelo=0 and Utelo=64. We can also examine the effects of knocking out different genes for checkpoint proteins and nucleases. For example, to simulate a cdc13-1 rad9Δ strain, we set the initial amount of Rad9I=0 (see table 4). Our model predicts that cdc13-1 cells are unable to divide, although occasionally they may divide once before cell-cycle arrest occurs (data not shown). This is due to the uncapped telomeres being bound by Rad17, which then initiates a series of events leading to cell-cycle arrest. Both Rad9 and Rad24 are required for cell-cycle arrest to occur. Our model predicts that cdc13-1 exo1Δ double mutants enter cell-cycle arrest after a few divisions as ssDNA can still be produced through the activity of ExoX and so Rad9 is activated, and in turn, the other checkpoint proteins are activated. The later activation of the checkpoint response is due to the slower generation of ssDNA without the activity of Exo1. Figure 6 shows data from Zubko et al. (2004) for various yeast strains cultured at the restrictive temperature. The cdc13-1 mutants are unable to divide at the restrictive temperature, but knocking out Rad24, Rad9 or Exo1 restores cell division, hence microcolonies form. Both cdc13-1 rad9Δ and cdc13-1 exo1Δ double mutants at the restrictive temperature have very small microcolonies, although we know that knocking out Rad9 would prevent cell-cycle arrest (figure 6). The cells for cdc13-1 exo1Δ strain are larger than in the cdc13-1 rad9Δ strain, which indicates that cells initially arrest and then eventually escape this arrest and divide a few times before cell death takes place. The cdc13-1 strains accumulate large amounts of ssDNA, since nuclease activity is intact and this is probably the cause of cell death. If we compare the model output with experimental data, we see good agreement apart from the cdc13-1 rad9Δ mutant. Our model predicts that cdc13-1 rad9Δ cells keep dividing, but experimental data show that these cells only form very small microcolonies (figure 6). This is because we had not included the possibility of cell death in the model. To take cell death into account, we modified the model (see §3.2.3). Our model also does not fit the observation of cell-cycle arrest followed by escape for the cdc13-1 exo1Δ strain as we have not included this mechanism in the model. We outline a possible extension to include this escape mechanism in §4.

Figure 6

Growth of wild-type and cdc13-1 mutant strains (reproduced with permission from Zubko et al. (2004); Copyright Genetics Society of America). Yeast strains were released from G1 arrest and allowed to form microcolonies for 15 h at 36°C (restrictive temperature) before being photographed at 200× magnification. Cell numbers within microcolonies were estimated from the photographs shown and are given, along with their standard deviations, for each strain: (a) cdc13-1 RAD+(2); (b) cdc13-1 exo1 (8±2); (c) cdc13-1 rad9 (18±5); (d) cdc13-1 rad9 ex01 (54±17); (e) cdc13-1 rad24 (43±22); (f) cdc13-1 rad9 rad24 exo1 (113±34); (g) wild type.

View this table:
Table 4

Summary table of knockout simulations for cdc13-1 strains (model with cell death).

3.2.2 Generation of ssDNA

Zubko et al. (2004) showed that large amounts of ssDNA are produced in cdc13-1 mutants at the restrictive temperature, but that deleting Exo1 considerably reduced these levels. Deleting Rad24 also reduced the amount of ssDNA, but deleting Rad9 had the opposite effect (see Zubko et al. 2004). Therefore, we plotted the level of ssDNA at different time points for each of the different mutant strains and compared our results to Zubko et al.'s (2004) data.

Our models predict that initially large amounts of ssDNA are generated in cdc13-1 strains (figure 7), but this stabilizes when cells enter G2/M arrest, as nuclease activity is inhibited by Rad9 and Rad53. Large amounts of ssDNA are generated in cdc13-1 rad9Δ double mutants, since Rad9 is required to inhibit nuclease activity. In contrast, only small amounts of ssDNA are generated in cdc13-1 rad24Δ double mutants. This is due to the fact that in our model Rad24 is required to load Rad17 onto uncapped telomeres and so activate ExoX. Our model predicts that much lower levels of ssDNA are generated in cdc13-1 exo1Δ double mutants than in the cdc13-1 strain. This low level of ssDNA is due to the residual activity of ExoX, which has not been completely inhibited by Rad9. The generation of ssDNA in the cdc13-1 rad24Δ strain is due to the Rad24-independent activity of Exo1 and the difference between the amount of ssDNA in the cdc13-1 strain and the cdc13-1 exo1Δ strain reflects the total activity of Exo1. Since experimental data show that this difference is larger than the amount of ssDNA generated in the cdc13-1 rad24Δ strain (Zubko et al. 2004), we suggest that there is also some additional Rad24-independent activation of Exo1 and so have included this in the model. Note that our model does not distinguish between different regions of the telomeres and that the model output represents the total amount of ssDNA for all telomeres. To examine the effects of nuclease activity at different regions of the telomere, we would need to modify the model to include a separate species for each telomere and separate reactions for different parts of the telomere.

Figure 7

Model predictions for the amount of ssDNA generated in wild-type and cdc13-1 mutant strains. The output is for one simulation for each strain over a period of 250 min.

We also looked at the variability of the amount of ssDNA produced over repeated simulations. We found that the variation from run to run was very small and even with 10 simulations, the upper and the lower 95% confidence interval for the mean was always within 2% of the mean. Therefore, it was not necessary to carry out many repeated simulations.

3.2.3 Extending the model to include cell death

When cdc13-1 rad9Δ cells are cultured at the restrictive temperature, it has been found that despite the lack of Rad9, which is required to activate the checkpoint pathway, these cells only divide for a few hours and then die (Zubko et al. 2004). The death of cells is due to the large amounts of ssDNA generated in these mutants. Therefore, we modified our model to allow for the possibility of cell death if the total ssDNA for the cell reached a critical threshold of 20 kb. There is currently no data on the amount of ssDNA required to trigger cell death, but it has been observed that there is ca 15 kb ssDNA present in cdc13-1 exo1Δ strains (Zubko et al. 2004). Hence, assumed that the threshold must be greater than this and initially chose 20 kb. The results for this model are in figure 8 and the model now predicts that the cdc13-1 rad9Δ cells do not divide since the critical threshold of ssDNA for cell death is reached before it is time to divide (figure 9a). Therefore, we increased the critical threshold for cell death to 120 kb and this allows the cdc13-1 rad9Δ cells to divide two to three times before cell death (figure 9b). These results agree much better with the experimental data, since the cdc13-1 rad9Δ exo1Δ cells are able to form larger microcolonies than both cdc13-1 rad9Δ cells and cdc13-1 exo1Δ cells (figure 6).

Figure 8

Model predictions for the number of divisions obtained by wild-type and cdc13-1 mutant strains if a critical threshold of 20 kb ssDNA triggers cell death. The output is for one simulation for each strain over a period of 12 h.

Figure 9

Model predictions for (a) the amount of ssDNA per cell and (b) the number of cell divisions if a critical threshold of 120 kb ssDNA triggers cell death. The output is for one simulation for each strain over a period of 12 h.

However, with this increased threshold, the model predicts that cdc13-1 rad9Δ exo1Δ cells are able to keep dividing as in the absence of Exo1, the rate of single strand production is lower and since these cells are dividing, ssDNA can be removed during S phase. This accounts for the zig-zag line in figure 9a. Since the rate of ssDNA removal depends on the level of ssDNA, there comes a point at which there is a balance between the rate of generation and removal of DNA. For the cdc13-1 rad9Δ exo1Δ cells, this balance occurs somewhere between 60 and 80 kb. This suggests that the threshold for cell death might not be as high as 120 kb, but somewhere in the region of 60–80 kb. Therefore, we did further simulations with a critical threshold of 70 kb. In this case, the model predicts that cdc13-1 rad9Δ exo1Δ cells divide 6–7 times, cdc13-1 exo1Δ cells divide 3–5 times, cdc13-1 rad9Δ cells divide just once and cdc13-1 cells do not divide at all (data not shown).

Simulations were also carried out for rad53Δ and mec1Δ mutant strains. We found lower rates of ssDNA production in both cdc13-1 rad53Δ exo1Δ and cdc13-1 mec1Δexo1Δ strains than in cdc13-exo1Δstrains, which agrees with experimental data (Jia et al. 2004) This ssDNA was the result of ExoX activity, which was not completely inhibited by Rad9. Our models predicted that the rate of ssDNA production was less than 1 kb h−1 on an average for both the strains (data not shown).

Table 4 summarizes all the knockout experiments performed and compares the model predictions with the experimental data. Where agreement is not close, a suggestion for modifying the model has been made.

3.3 Deterministic versus stochastic simulation

We also used a deterministic simulator to run the model to examine whether any parts of the model could be successfully modelled using an ordinary-differential-equation (ODE) framework. We found that the generation of ssDNA could be modelled in this way, and the output for the amount of ssDNA is very similar using either a deterministic or stochastic simulator. This is because the variation in the amount of ssDNA for a number of stochastic simulations was small and typically large numbers were involved. However, information was lost for other parts of the model when using a deterministic simulator. For example, in the wild-type model, the deterministic output for the number of capped and uncapped telomeres gave a steady state of 63.85 and 0.15, respectively, which reflects the average number of capped and uncapped telomeres at any point in time, whereas output from the stochastic simulator show the exact numbers of capped and uncapped telomeres at any time point. Similarly, the number of divisions obtained by cdc13-1 mutant strains in a stochastic simulation was generally low integers and showed some variability, whereas the deterministic output could only show the average value. The advantage of using deterministic models is the greater speed of simulation and so it is worth considering the usage of a hybrid approach, where parts of the model are simulated using a deterministic simulator and other parts (e.g. low-copy number species) are simulated with a stochastic simulator. There are also hybrid simulators being developed, which combine exact stochastic with approximate stochastic algorithms. Further details of hybrid simulation can be found in Kiehl et al. (2004) and Wilkinson (2006).

4. Discussion

We have developed a mathematical model of the yeast checkpoint response, which has been successfully used to test ideas about the checkpoint response to telomere uncapping in budding yeast. The model was modified and used to predict the behaviour of cdc13-1 cells at the restrictive temperature, and in addition, we carried out simulations, where one or more genes were deleted. We compared simulation output to experimental data of the amount of ssDNA generated over time and the number of cell divisions obtained for the different mutant strains. Generally, we found good agreement with our model predictions and the schema in figure 1. We found that it was necessary to include cell death in the model and that the model output depends on the assumption made about the critical threshold of ssDNA for cell death. This has been useful to discover and suggests that further experimental work is needed to determine this threshold. To explain the experimental observation that cdc13-1 exo1Δ initially stops dividing and then escapes from cell-cycle arrest, it will be necessary to extend the model to include the mechanism of adaptation (discussed later).

The models have been kept as simple as possible while retaining enough biological realism and provide the starting point for further investigation. They were encoded in SBML, hence they can easily be extended at any point in the network. For example, it is a fairly easy task to add in further proteins that are involved in the network such as the MRX complex (Mre11, Rad50 and Xrs2), which may be involved in processing telomeric ends and in telomere capping (Larrivee et al. 2004; Foster et al. 2006). Another protein complex, important at telomeres, is the Yku70/Yku80 heterodimer, which has a similar role to Cdc13 in binding telomere ends and protecting them from repair and checkpoint pathways (Maringele & Lydall 2002). An important kinase which we have not yet included in the model is the Tel1 kinase. Takata et al. (2004) show that Tel1 and Mec1 are recruited reciprocally to telomeres during the cell cycle. Tel1 is recruited to telomeres with a repressive structure and is needed to prevent telomeres from fusing through non-homologous end joining in collaboration with telomerase. Mec1 associates preferentially with shortened telomeres during replication. Our model could be modified to investigate this process.

It might also be desirable to add in further details about some of the proteins that are already in the model. For example, Rad9 functions as an 859 kDa complex in undamaged cells, but undergoes a conformational change in the presence of DNA damage (Gilbert et al. 2001). This change involves loss of mass and hyperphosphorylation and requires the essential chaperones, Ssa1 and Ssa2 (Gilbert et al. 2003). Another example is Mec1, which not only is recruited to damaged DNA, but also resides at the telomeres even in the absence of damage (Takata et al. 2004). This might affect the kinetics of RPA recruitment to ssDNA and therefore it requires further investigation. Interestingly, a recent study shows that the cell-cycle kinase Cdk1 is required for the recruitment of ssDNA-binding complex, RPA and that it is also involved in degradation, which occurs at double-stranded DNA breaks (Ira et al. 2004). There is no direct evidence that it is involved in the degradation of ssDNA at telomeres, but it would be interesting to see how adding the requirement for Cdk1 to allow RPA-binding affects the model predictions.

Further details of the telomere end protection pathway could also be incorporated into the model. For example, a recent study by Verdun et al. (2005) found that telomeres are unprotected during the G2 phase of the cell cycle and that a localized DNA damage response at telomeres after replication is essential for recruiting the processing machinery that promotes formation of a chromosome end protection complex. Since our model incorporates the cell cycle, it would be fairly straightforward to add this extra detail to the model.

A less straightforward, but by no means impossible, task is to replace the pool of telomeres with individual telomeres to allow the investigation of telomere dynamics. The easiest way to do this would be to have an array or a set of telomere species in the model. Currently, SBML does not support models with arrays or sets, but proposals have been put forward to incorporate this feature into SBML level 3 (www.sbml.org). The model could then be used to investigate the action of telomerase and the possibility that either abrupt shortening or lengthening may cause cell-cycle arrest (Ijpma & Greider 2003; Viscardi et al. 2003). To examine the effects of nuclease activity at different regions of the telomere (Zubko et al. 2004), we would need additionally to add separate reactions for different parts of the telomere and keep track of the amount of ssDNA acquired in the different regions. This would be very interesting to do and the model could then be used to test the hypotheses that different nucleases are important at different regions of the telomeres (Zubko et al. 2004).

We could also modify our model so that it is possible to carry out simulations, which involve changing the temperature of cdc13-1 mutant strains. A simple way to do this is to add two variables to the model to represent the permissive and restrictive temperature, e.g. Temp23 and Temp36. At the permissive temperature, Temp23 and Temp36 are set to 1 and 0, respectively, and at the restrictive temperature, the settings are reversed. The addition of Temp23 to the kinetic law of the capping reaction would ensure that this reaction could only proceed at the permissive temperature (if Temp23=0, then the rate of the reaction would also be 0 and so unable to take place). The temperature could be switched at certain time points during the course of the simulation to mimic the experimental procedure of moving cells from the permissive temperature to the restrictive temperature and vice versa.

Experimental data shows that the function of cdc13-1 is totally disrupted at 36°C (Garvik et al. 1995). As the temperature is lowered, cdc13-1 function is gradually restored so that there is a positive correlation between function and decreasing temperature (Zubko et al. 2004). We could model this scenario by having a species in the model to represent the temperature, e.g. Temp, and to choose a kinetic law so that the rate of the capping reaction depends on the value of Temp. The value of Temp could be changed at different time points during the simulation to mimic the experimental procedure of switching temperatures.

The cells for cdc13-1 exo1Δ strain are larger than in the cdc13-1 rad9Δ strain, which indicates that cells initially arrest and then eventually escape this arrest and divide an average of three times before cell death takes place (figure 6; Zubko et al. 2004). Our model predicts that the cdc13-1 exo1Δ strain divide four or five times and then arrest. Therefore, it is desirable to add the possibility of escape from cell arrest into the model. It is known that escape (also known as adaptation) from cell arrest occurs even though the damage and the signal for damage persists and that the time for escape is ca 8–15 h (Toczyski et al. 1997; Gardner et al. 1999; Pellicioli et al. 2001). Studies have shown that escape requires the proteins Cdc5 and casein kinase II (Toczyski et al. 1997) and that Rad53 kinase activity and Chk1 phosphorylation largely disappear at the time that cells escape arrest, provided functional Cdc5 is present (Pellicioli et al. 2001). Cdc5 is a target of Rad53, and phosphorylation of Cdc5 is required for the completion of anaphase (Sanchez et al. 1999). Pellicioli et al. (2001) suggest that Cdc5 acts in a feedback loop to turn off the checkpoint kinase cascade, but at present we do not know the exact mechanism of how this is achieved. Trying to model this scenario would help to clarify the details of this pathway.

The models in this paper have followed the fate of one individual cell. It would be interesting to extend the models to follow an entire colony. This would be much more computer intensive to do, but should become easier in the future as computing power continues to increase. The addition of arrays or sets into SBML would also make it much easier to code population-based models.

Finally, the models developed in this paper were motivated by hypotheses put forward as a result of experimental work. Our models do not preclude the possibility that other hypotheses may also be able to explain the experimental results. Therefore, we strongly encourage the interested reader to either make their own modifications to the models or to send us their suggestions so that other ideas can be tested.

Acknowledgments

This work was funded by BBSRC, MRC, EPSRC, DTI and Unilever plc. D.L. is also funded by the Wellcome Trust. We thank three anonymous referees for their very helpful suggestions.

Footnotes

    • Received June 27, 2006.
    • Accepted July 14, 2006.

References

View Abstract