Model
The model consists of population of replicating and dying cells. The population is limited to 400 replicating cells: when the population size drops below 400 a random cell, that has divided more than 24 hours (set to be cell doubling time) prior to present event, is picked and is allowed to divide. In the absence of genotoxic stress, cells are dying stochastically (due to DNA damage unrelated reasons) with a constant rate β
_{0}=0.02 such that on average each cell doubles every 24 hours. To easier relate our model to the experimental data we choose to report our results in units of “cell doublings” with one cell doubling being equivalent to 24 hours in our model. The population starts with 400 replicating cells having initial telomere length of 15000 bp. Cells are removed from the pool of replicating cells if either telomere length is below critical length (2000 bp) or the p53 mediated mechanisms are activated. In the following sections we present the biological basis for our model, followed by the detailed description of the model implementation.
Biological basis of the model
The main components of the model are: i)DNA damage, D, ii) accumulation of mutations, M, iii) cell cycle arrest, apoptosis and other p53 mediated responses, p 53 and iv) telomere shortening, T (see Figure 1).
i) The level of DNA damage, D: is assigned from a normal distribution with a given average μ _{
D
} and a standard deviation σ _{
D
}. The variation in D might arise from multiple sources: fluctuations in concentrations of DNA repair enzymes, variations in metabolic load (and subsequent Reactive Oxygen Species (ROS) production) in single cells, etc. The particular choice of the shape of the distribution (e.g. it has been shown that celltocell variation in gene expression sometimes follows lognormal distribution [10, 11]) does not qualitative change the main outcomes of the model.
As the timescale relevant for our simulation is of the order of cell division (characteristic timescale for both telomere decay and mutation rate) we assign a new value of D every time cells divide. In fact any fluctuations much faster than cell doubling time will be averaged out and result in a “homogeneous population” where each cell experiences the same damage seen at the timescale of cell doubling. The other limit, when fluctuations in D are much slower than doubling time, will again result in somewhat “homogeneous population” with several different groups of cells. Thus the most interesting regime is when D changes on timescale of cell doubling. Note that while we are assigning the damage from Gaussian distribution, the resulting distribution of damage in the simulated population of replicating cells can be different from Gaussian (e.g. damage can not be negative). In the following we will denote the damage averaged over cells and time as 〈D〉.
ii) Mutations, M spread and accumulate in the population as mutated cells replicate. Cells with many mutations have higher chance to originate tumor cells [12] and can be thought of as tumor progenitors.
The genotoxic stress, e.g. oxidative stress, replication of DNA fragile sites, UV or gamma radiation, etc. results in DNA damage (D, Figure 1). Damaged DNA recruits DNA repair machinery by e.g. activating ATM signaling cascade. In most of cases DNA repair enzymes remove the damage [13–15], however the repair is not perfect and often mutations occur as a result of damage and repair cycles. Thus higher genotoxic stress leads to more mutations [5]. We model this dependence by setting the rate of probability for a mutation to occur,
{r}_{M}=\alpha \frac{{D}^{2}}{{D}^{2}+1}.
(1)
The typical mutation rate is estimated to be 10^{11} (somatic stem cells)  10^{9} (typical for proliferated cells) per basepair per cell division, which amounts to 0.011 per human cell per cell division [12]. As we are simulating observations made in proliferated cells, we set α=1 per cell division. We choose mutation rate of r
_{
M
}=0.2 per cell per cell division to represent “typical” mutation rate under “physiological” damage (thus with α=1 the range of “physiological” damage is 〈D〉∼0.5). The main results will be qualitatively the same if the mutation rate is increased or decreased 5 fold.
iii) p 53mediated responses limit the mutation spread by rapidly (compared to the telomere attrition) eliminating stressed cells. Cell survival under DNA damage have sigmoid doseresponse curve [16], we have modeled this observations by setting probability for p53 mediated cell death to be a sigmoid curve
{r}_{\mathit{\text{apopt}}}=\beta \frac{{D}^{2}}{{D}^{2}+1}
(2)
Observe that the functional form in mutational and apoptotic probabilities are set to be the same as this allows most efficient elimination of mutated cells. Parameter β=0.1 is chosen such that just a small fraction of cells (0.1%) undergoes apoptosis at low levels of D∼0.2.
The results of the model do not depend on the choice of the functional forms for probabilities to mutate or undergo apoptosis. (See Additional file 1: Figure S1.)
iv) Telomeres consist of (TTAGGG) repeats which form a protective cap at the end of eukaryotic chromosomes. During cells division, the 3’ end of linear DNA can not be fully replicated and thus telomeres become shorter. Cells with short telomeres (2000 bp) lose the ability to replicate and at this stage they either exist in a nondividing state or undergo programmed cell death. Interestingly, the rate of loss of telomeric DNA is not constant but appears to depend on the length of telomere [17, 18] and the level of oxidative stress [9]. In our model we capture these observation by setting the telomere decay to be proportional to DNA damage, D[8] and the length of telomere, T[17, 18].
\frac{\mathit{\text{dT}}}{\mathit{\text{dt}}}=\mathrm{\gamma TD}
(3)
Thus we do not explicitly model the mechanism of how telomeric damage leads to telomere decay. This has been carefully addressed in the model by Proctor et al. [8] and is beyond the scope of our work. Instead we phenomenologically describe the observed correlation between the rate of telomere decay and DNA damage in the cell and assume that the cellular damage is independent of telomere length. While we model one telomere per cell, in reality there 92 telomeres per cell. Telomere lengths follow a skewed, lognormallike, distribution and it is believed that the the replicative senescence is dictated by the shortest telomere [19]. Furthermore, we model telomere decay as a continuous process while in living cells the decrease in telomere length is related to the replication and happens at cell division. Replacing continuous update of telomere lengths with a discrete update leads to same qualitative results.
Initial telomere length in human fibroblasts was estimated to be 15000 bp and the rate of decay is about 100 bp per division [7]. The exact values of the initial and critical telomere lengths do not affect the qualitative results of the model. Parameter γ=1.5×10^{3} is given by the requirement of decay of 50100 bp per cell division when DNA damage, D is low (D=0.2). This description of telomere dynamics is inspired by the model by Proctor et al. [8]. As parameters α,β and γ are constrained by experimental data, the only free parameters of interest are average DNA damage, 〈D〉 and how the damage differs from celltocell, σ
_{
D
}.
Model execution
The code executing the model is programed in C++ and completes within minutes on a standard PC. For each cell we keep track of the following attributes:
D DNA damage, assigned at cell division from gaussian distribution with μ
_{
D
} and σ
_{
D
}
t
_{
birth
} birth time, set for daughter cells after every division
τ
_{
surv
} survival time, assigned at cell division. To arrive to a damageinduced rate of cell death given by Eq. 2 and account for stochastic damageunrelated death with rate β
_{0}, τ
_{
surv
} is drawn from exponential distribution {e}^{{\beta}_{0}\beta {D}^{2}/({D}^{2}+1)(t{t}_{\mathit{\text{birth}}})}.
T Telomere length, updated every time step.
n
_{
mut
} mutation counter, updated every time step.
The time evolution of the model is as follows:
0) At time t=0, for each cell among 400 cells we initialize T=15000, n
_{
mut
}=0, t
_{
birth
}=0 and assign D and τ
_{
surv
} as described above.
1) At every time step advance in time with d t=0.1h and in each cell

Update T, according to Eq.3.

Increase n _{
mut
} by one with probability given by Eq. 1

If T≤2000 or t _{
current
}t _{
birth
}≥τ _{
surv
}, i.e. if cell turns senescent or undergoes apoptosis
* remove cell from the population of dividing cells.
* Divide a cell chosen randomly among those with t
_{
current
}t
_{
birth
}≥24.
* Daughter cells inherit T,n
_{
mut
} and get assigned new D and τ
_{
surv
} as described above. For each of the daughter cells set t
_{
birth
}=t
_{
current
}.
2) Repeat advancing in time as described in 0) until there are no replicating cells left in the population.