Cell cycle progression

Joanna Tyrcha

doi:10.1016/j.crvi.2003.05.002

Résumé

In this paper we consider cell cycle models for which the transition operator for the evolution of birth mass density is a simple, linear dynamical system with a stochastic perturbation. The convolution model for a birth mass distribution is presented. Density functions of birth mass and tail probabilities in n-th generation are calculated by a saddle-point approximation method. With these probabilities, representing the probability of exceeding an acceptable mass value, we have more control over pathological growth. A computer simulation is presented for cell proliferation in the age-dependent cell cycle model. The simulation takes into account the fact that the age-dependent model with a linear growth is a simple linear dynamical system with an additive stochastic perturbation. The simulated data as well as the experimental data (generation times for mouse L) are fitted by the proposed convolution model.

1 Introduction

The cell cycle is normally defined as the period from cell birth to its division into two daughter cells. The most significant aspect of cell growth is the synthesis of new cellular material, such as proteins, and the initiation of DNA synthesis is one of the first observable cell cycle landmarks. There are cells that are characterized by uncontrolled growth. It is of greatest interest to find answer to the question what causes the cell to begin initiation of DNA synthesis, since this problem is closely connected with explaining the cause of pathological growth, such as in cancer.

Many mathematical models were proposed which explain the asymptotic stability of the cell cycle. By the stability concept we mean the stability of the line: mother–daughter–granddaughter–great-grand-daughter, etc. We derive a recursion relation for the distribution of cell mass at birth in a sample of cells in the generation n given the distribution in the generation n−1. We prove the existence of a unique, asymptotically stable, limiting mass distribution. An alternative convention suggested, among others, by Harris [1] and Jagers [2] is to choose samples from all cells present in the culture at a particular instant of time, and investigate stabilization of the mass distribution in the population as time passes.

The cell cycle is one of many examples of systems that have a mixture of deterministic and probabilistic dynamics. Many other examples may be found in [3]. In this paper we develop a general modelling framework for such systems.

First, we shortly review previous cell cycle models and the mathematical tools used to investigate the asymptotic stability of the mass distribution at birth. We concentrate on the models that represent linear dynamical systems with a stochastic perturbation. In Section 3, we present the convolution model for estimating distributions of cell mass at birth. In Section 4, we describe our simulation procedure. In Section 5, we present the maximum likelihood method for estimate parameters in our convolution model and finally show how the theoretical density functions of cell mass birth fit both our simulated and experimental data.

2 Background

The cell mass at birth (x) may be considered as a random variable. If we denote by f the density function of the mass distribution at birth for a mother cell then Pf denotes the density function of the mass distribution at birth for daughter cells. In general the operator P has the form

Pf (x) = \underset{0}{\int^{\infty}} K (x, y) f (y) d y

(1)

where K is to be thought of as the conditional probability that, given a mass y in the previous generation, we will have a mass x in the current generation. K is a stochastic Markov kernel, which means:

\underset{0}{\int^{\infty}} K (x, r) d x = 1 for r ⩾ 0 and

K (x, r) ⩾ 0 for x ⩾ 0, r ⩾ 0

The particular form of the kernel K(x,y) depends on the assumptions concerning the mass growth and the probability of division. Depending on these assumptions we can distinguish, among others, the following models of the cell cycle: the mitogen model [4], the tandem model [5], the probabilistic model [6,7], and the generalized mitogen model [8]. In all these models, the authors studied asymptotic properties of the sequence of iterates {Pⁿ} (Pⁿ denotes the mass density function in the n-th generation).

In the paper by Lasota, Mackey and Tyrcha [9], we developed a general modelling framework for the treatment of the statistical dynamics of systems in which easily identifiable events occur at irregular times. The desired recurrence relation between successive mass levels at event occurrence is given by:

x_{n + 1} = λ^{- 1} (Q^{- 1} (Q (x_{n}) + τ_{n})) for n = 0, 1, \dots

(2)

where, by assumption, x_n and τ_n are independent stochastic variables. We assume that Q and λ satisfy the following conditions:

– the functions $Q : R_{+} \to R_{+}$ and $λ : R_{+} \to R_{+}$ are non-decreasing and absolutely continuous on each subinterval [0,c] on the half-line R₊;
– Q(0)=λ(0)=0 and $\lim_{n \to \infty} Q (x) = \lim_{x \to \infty} λ (x) = \infty$ .

In interpreting the cell division cycle within the context of the general model presented in [9], we associate the occurrence of an event with a triggering of the process, which ultimately leads to mitosis and cytokinesis.

The class of models proposed by Lasota and Mackey [4], Tyson and Hannsgen [5] and Tyrcha [8] satisfies all of the conditions of the general model described above and in [9]. Furthermore, the quantities of mass in consecutive generations of newly born cells satisfy the recurrence relation (2) with τ_n having a survival function e^−x. Setting z_n=Q[λ(x_n)], the transition operator for the evolution of mass density (2) is given by:

z_{n + 1} = Q [λ^{- 1} (Q^{- 1} (z_{n}))] + τ_{n}

(3)

As another example, extensions proposed by Tyson and Hannsgen [7] and Hannsgen et al. [6] of the well-known cell cycle models of Smith and Martin [10] and Shields [11] also fall within the general framework of this paper. In this situation, we assume that the cell goes through phases A and B, and that the length T_B of phase B is constant. The end of the B phase marks cell division. The length T_A of phase A is considered to be a random variable with a density function h so as:

Prob (T_{A} ⩾ t) = \underset{t}{\int^{\infty}} h (z) d z for t ⩾ 0

(4)

The mass is produced with dynamics described by equation dx/dt=g(x), g>0, is assumed to not affect T_A, and divides equally between daughter cells at division. Thus, by assumption, T_A and x(0) (x(0) is a cell mass at time t=0) are independent. So, the dynamical system (3) reduces to [12]:

z_{n + 1} = \frac{1}{2} z_{n} + \frac{1}{2} τ_{n}

(5)

where z₀ is the cell mass at birth in 0-th generation and τ_n (n=0,1,…) is a stochastic mass increment during the time T_{A_n}+T_B. This is a simple, linear dynamical system with a stochastic perturbation.

3 Convolution model

We considered before [8,12] the asymptotic stability of the cell mass at birth for both mass-dependent and age-dependent models of the cell cycle. In this paper, we describe a cell mass distribution at birth in nth generation if we assume a cell mass distribution at birth in 0-th generation and distribution of stochastic perturbations (Eqs. (3) and (5)).

Let x₀, τ₀, τ₁, ‘K’ be a sequence of independent stochastic variables where x₀ is a Gaussian stochastic variable and the τ_n are exponentially distributed. Our model for the first cell generation is a truncated distribution of the convolution of the Gaussian and the exponential distributions. Such a model is a result of discussion with a group of biologists lead by Prof. Anders Zetterberg (Karolinska Institute, Stockholm, Sweden). From experimental investigations, it follows the assumption about initial cell mass distribution at birth being Gaussian (see also [13]). Starting from the initial mass, a cell is growing up but, even if it grows up according to a deterministic equation, one can observe some stochastic perturbations. A stochastic mass increment is assumed to be exponentially distributed. We are interested in the truncated distribution, because we would like to eliminate cells with an extremely high mass.

The distribution of the convolution of the Gaussian and the exponential distributions depend on three parameters: the mean value (equal to the standard deviation) 1/a of the exponential, where a is the intensity parameter, and the mean μ and the standard deviation σ of the Gaussian component. Its density can be expressed as:

f (x) = a \exp \{\frac{a^{2} σ^{2}}{2} - a (x - μ)\} Φ (\frac{x - μ - σ^{2} a}{σ})

(6)

where Φ denotes the standard Gaussian cumulative distribution function.

When cell masses above a maximum value T are not taken into consideration, it means that distribution in Eq. (6) is truncated in T. The truncated density function f_T(x) looks as follows:

f_{T} (x) = {\begin{matrix} \frac{f (x)}{F (T)}, & 0 ⩽ x ⩽ T \\ 0, & otherwise \end{matrix}

(7)

where

F (T) = Φ (\frac{T - μ}{σ}) - \exp \{\frac{a^{2} σ^{2}}{2} - a (T - μ)\} Φ (\frac{T - μ - σ^{2} a}{σ})

is the cumulative distribution function corresponding to f.

If we assume that Eq. (7) holds, then our model for the n-th cell generation is the truncated convolution of n copies of the distribution for one generation. The convolution of f_T with itself cannot be explicitly given. Usually calculations for convolutions are made via Fourier transforms and Fourier inversion, since a convolution of distributions corresponds to multiplication of their Fourier transforms. However, the truncated distribution (7) for the n-th cell generation has no explicit Fourier transform.

Instead, we have calculated the distribution for the n-th cell generation by a saddle-point approximation [14–16]. This is a very accurate distribution approximation in many situations. The saddle-point approximation method uses the Laplace transform (the moment generating function), and not the Fourier transform. The Laplace transform ψ(θ) of (6) can be explicitly calculated:

ψ (θ) = \frac{1}{F (T) (1 - θ / a)} [b_{1} (θ) - b_{2} (θ)]

where

b_{1} (θ) = \exp \{μ θ + \frac{σ^{2} θ^{2}}{2}\} Φ (\frac{T - μ - σ^{2} θ}{σ})

and

b_{2} (θ) = \exp \{\frac{a^{2} σ^{2}}{2} + a μ + T (θ - a)\} Φ (\frac{T - μ - σ^{2} a}{σ})

Let X₁,…,X_n be independent and identically distributed random variables, copies of a random variable X with density f_T(x) and Laplace transform ψ(θ). Consider:

S_{n} = X_{1} + \dots + X_{n}

The density function f_n of S_n is the n-fold convolution of f_T:

f_{n} (s) = f_{T}^{n^{*}} (s)

As a technical trick (exponential tilting), we associate with f_n(s) the exponential family:

f_{n} (s; θ) = e^{θ s} f_{n} (s) / ψ {(θ)}^{n} = \exp \{θ s - K (θ)\} f_{n} (s)

(8)

where K(θ)=nlogψ(θ) is the cumulant transform (the cumulant generating function) for S_n. We have the formulas

μ (θ) = \frac{\partial K}{\partial θ}, σ^{2} (θ) = \frac{\partial^{2} K}{\partial θ^{2}}

C_{m} (θ) = \frac{K^{(m)} (θ)}{σ^{m} (θ)} (m > 2)

where μ(θ) and σ²(θ) denote respectively the expected value and the variance for the exponential family f_n(s;θ), and K^(m) and C_m are raw and normalized cumulants, respectively. These formulae are well known from the general theory of exponential families, see, for example, [14,17]. In particular, note that μ=μ(0) and σ=σ(0). More explicitly, in our case:

μ (θ) = \frac{1}{a - θ} + \frac{(μ + θ σ^{2}) b_{1} (θ) - {Tb}_{2} (θ) - b_{3} (θ)}{b_{1} (θ) - b_{2} (θ)}

(9)

where the functions b₁ and b₂ are given above and

b_{3} (θ) = \frac{σ}{\sqrt{2 π}} e^{θ T - {(T - μ)}^{2} / (2 σ^{2})}

The second derivative of the cumulant transform K has the form:

σ^{2} (θ) = \frac{1}{{(a - θ)}^{2}} + [σ^{2} b_{1} (θ) + (μ + θ σ^{2}) \{(μ + θ σ^{2}) b_{1} (θ) - b_{3} (θ)\} - T^{2} b_{2} (θ) - {Tb}_{3} (θ)] \cdot {(b_{1} (θ) - b_{2} (θ))}^{- 1} - \frac{{(μ + θ σ^{2}) b_{1} (θ) - {Tb}_{2} (θ) - b_{3} (θ)}^{2}}{{b_{1} (θ) - b_{2} (θ)}^{2}}

We require a good approximation to the tail probability Prob(S_n>s) because it is interesting to be able to calculate the probability that a cell mass exceeds some level. However, we start by giving an approximation to the density f_n of S_n, which is analogous, but slightly simpler. For any such s, the corresponding root

\hat{θ} = \hat{θ} (s)

of the equation K′(θ)=s (the saddle point) has a crucial role. Suppose now that we wish to calculate a good approximation to f_n(s)=f_n(s;0) for a given value of s. We write:

f_{n} (s) = \exp \{- θ s + K (θ)\} f_{n} (s; θ)

and then choose

θ = \hat{θ} (s)

, that is the ML estimator in the exponential family (8). The approximation

{\hat{f}}_{n} (s) = \frac{ψ {(\hat{θ})}^{n} \exp (- \hat{θ} s)}{σ (\hat{θ}) \sqrt{2 π n}}

is the saddle-point approximation to the density function f_n of S_n. A refined version is:

{\hat{f}}_{n} (s) = \frac{ψ {(\hat{θ})}^{n} \exp (- \hat{θ} s)}{σ (\hat{θ}) \sqrt{2 π n}} \{1 + \frac{1}{n} [\frac{C_{4} (\hat{θ})}{8} - \frac{5}{24} C_{3} {(\hat{θ})}^{2}]\}

see [14] (formula (2.2.4)). Our model of the cell birth mass in n-th generation is the truncated convolution of n copies of the distribution given by (7), so the density function

{\hat{f}}_{n}^{*} (s)

for the truncated s is:

{\hat{f}}_{n}^{*} (s) = {\begin{matrix} \frac{{\hat{f}}_{n} (s)}{F (Tr)}, & 0 ⩽ s ⩽ Tr \\ 0, & otherwise \end{matrix}

where F(Tr) is the cumulative distribution function corresponding to

{\hat{f}}_{n} (s)

. Before presenting the formula for the cumulative distribution function F(s) or tail probabilities, we would like to make the following remark.

In many applications, the saddle-point equation μ(θ)=s cannot be solved analytically, although the solution $\hat{θ}$ exists. This is also the case here. The saddle-point method can still be applied by solving the equation numerically. Often one uses Newton–Raphson derivative-based methods to calculate the saddle-point. These methods in general will work well, since the function K(θ)−θs to be minimized is convex. The second derivative exists and is always positive. An initial approximation θ₀ is chosen and an iterative procedure is then used to calculate the saddle-point.

Let us turn then to the tail probabilities. We consider Prob(S_n⩾s), where s⩾μ(0). Then the solution $\hat{θ}$ of equation μ(θ)=s satisfies $\hat{θ} ⩾ 0$ . For $\hat{θ} ⩾ 0$ , we can calculate the tail probability (or risk) that the cell mass exceeds any specified level s according to the formula

Prob (S_{n} ⩾ s) = \frac{ψ {(\hat{θ})}^{n} \exp (- \hat{θ} s)}{\sqrt{n} | \hat{θ} | σ (\hat{θ})} B_{0} (η)

(10)

or the refined version

Prob (S_{n} ⩾ s) = \frac{ψ {(\hat{θ})}^{n} \exp (- \hat{θ} s)}{\sqrt{n} | \hat{θ} | σ (\hat{θ})} \{B_{0} (η) + \frac{sgn (\hat{θ})}{\sqrt{n}} \frac{C_{3} (\hat{θ})}{6} B_{3} (η)\}

where B₀(η) and B₃(η) are the Esscher functions:

B_{k} (η) = \frac{1}{2 π} \underset{- \infty}{\int^{\infty}} \exp (- \frac{t^{2}}{2}) \frac{{(i t)}^{k}}{1 + i t / λ} d t

see [14] (formula (2.2.6)). The Esscher functions can be expressed in terms of Φ by:

B_{0} (η) = η \exp (\frac{η^{2}}{2}) [1 - Φ (η)]

and

B_{3} (η) = - [η^{3} B_{0} (η) - (η^{3} - η) {(2 π)}^{- 1 / 2}]

Observe that the Edgeworth approximation methods provide reasonably good approximations in the centre of the density, but not in the tails, where the approximating density can even be negative. The saddle-point method gives good approximations farther out in the tail and quite good approximations to small-tail probabilities.

4 Data simulation

To see how our analytical formula for ${\hat{f}}_{n}^{*} (s)$ compares with a birth-mass histogram from a finite cell population, we have simulated an ideal cellular population by Monte Carlo methods. In particular, 550 cells were allowed to grow linearly and divide according to the assumptions in the age-dependent model. We chose the age-dependent model because it was shown [12] that together with a linear growth, this model is a simple linear dynamical system with an additive stochastic perturbation (see formula (5)). So, we assume that at division each mother cell divides equally between two daughter cells. The length T_A of phase A is a random variable exponentially distributed (formula (4)). At division, one daughter was removed from the culture, so that there were always 550 cells in each generation. From an initial arbitrary distribution of birth masses of the 550 cells, the culture was simulated for 20 generations. A histogram of 550 cell birth masses is shown in Fig. 1. The analysis presented in this paper assumes that the fundamentally discrete histogram of cell birth masses in a finite population can be approximated by a continuous probability density function ${\hat{f}}_{n}^{*} (s)$ .

5 Parameter estimation

5.1 Maximum likelihood method

We here discuss the estimation of model parameters from a sample of cell masses, specifically the sample illustrated in Fig. 1. We first assume that the model is given by (6), that is a three-parametric convolution of an exponential distribution with a Gaussian one. In the next section we allow occasional high masses, occurring with a probability δ, which is then an additional parameter to be estimated.

Let the sample be x₁, x₂, ‘Λ’, x_N. The likelihood function is:

L (a, μ, σ) = \prod_{i = 1}^{N} f (x_{i}, a, μ, σ)

where

f (x_{i}, a, μ, σ) = a \exp \{\frac{a^{2} σ^{2}}{2} - a (x_{i} - μ)\} Φ (\frac{x_{i} - μ - σ^{2} a}{σ})

(11)

The maximum likelihood estimates satisfy the likelihood equations:

\frac{\partial \log L}{\partial a} = 0, \frac{\partial \log L}{\partial μ} = 0, \frac{\partial \log L}{\partial σ} = 0

We have, taking natural logarithms,

\log L (a, μ, σ) = N \log a + N \frac{a^{2} σ^{2}}{2} - a \sum_{i = 1}^{N} x_{i} + Na μ + \sum_{i = 1}^{N} \log Φ (\frac{x_{i} - μ - σ^{2} a}{σ})

and we obtain the following equations: (12)

where

z_{i} = \frac{x_{i} - μ - σ^{2} a}{σ}

These equations, which can be simplified, are expressed above in the form E(·|x₁,…,x_n)=E(·), for a representation of the minimal sufficient statistic of the ‘complete’ exponential family we would have had if we could have observed the exponential and the Gaussian separately. It is impossible to solve these equations explicitly, so some iterative method is required. One possibility is to use the EM algorithm [18,19], which will iterate between the left- and right-hand sides of (12), but it will be quite slow. Therefore, we instead eliminated μ $(μ = \bar{x} - 1 / a)$ from the first equation and obtained the following profile-likelihood function for (a,σ):

\log L (a, σ) = N \log a + \frac{{Na}^{2} σ^{2}}{2} - Na \bar{x} + Na (\bar{x} - \frac{1}{a}) + \sum_{i = 1}^{N} \log Φ (\frac{x_{i} - (\bar{x} - 1 / a) - σ^{2} a}{σ})

We calculated the above log-likelihood function over a grid of different values of the parameters a and σ. The maximum point was found by making the grid successively finer. Based on the complete cell mass data shown in Fig. 1, the parameters were estimated to be: a=11.6, μ=0.6656, σ=0.0079. This estimated distribution did not fit the data well (solid line in Fig. 2).

Fig. 2
Histogram and two density functions.

5.2 Modified maximum-likelihood method

The following modified ML iterative procedure was found to be practical and adequate for dealing with the high cell masses in the large sample of data that was available (N=550). As the first step, model (6) is fitted by ML to all of data, as described above. In the fitted model, the (high) value t is determined such that only one observation should have been expected above t, according to (6). In our case, an explicit and close approximation to t is given by:

\exp (- \hat{a} (t - \hat{μ}) + {\hat{a}}^{2} {\hat{σ}}^{2} / 2) = 1 / N

because in the tail of (6) the Φ-factor will be close to 1. Now, all but one of the d (say) actual observations above t are assumed to represent what we call high losses, and consequently the probability δ for high losses can be estimated as

\hat{δ} = (d - 1) / N

. If exactly one observation were to be found above t, its expected value in the fitted distribution would be approximately

t + 1 / \hat{λ}

. Hence, as the second step, the d actual observations above t are replaced by a single observation at

t + 1 / \hat{λ}

, and the parameters of (6) are re-estimated.

This procedure is repeated until convergence.

With the actual data shown in Fig. 1, the parameters were estimated to be: a=15.0, μ=0.6673, σ=0.0089. Fig. 2 shows the fit of this distribution to the data truncated at 1.05.

6 Generation times

In the previous sections, we considered the cell mass at division for models of the cell cycle. Knowing the distribution of birth masses, we can also determine the distribution of generation times, i.e. the time from cell birth to division. In the particular case when $f_{n} = f_{*} (n = 0, 1, \dots)$ is a time-independent stationary sequence, the distribution of generation time has the same property.

We tried to verify our convolution model using experimental data from mouse L. We received the data from Prof. Michael C. Mackey (McGill University, Montreal, Canada). The data contain 89 pairs of generation times for respective daughter (sister) cells. Using the maximum-likelihood technique described in Section 5.1, we estimate parameters in the convolution model to be: a=3.2, μ=0.5178, σ=0.02. The theoretical density function fits the experimental data quite well (see Fig. 3). The same experimental data (mouse L) was used before to check the sister cell model described in [12]. The fitting of the data by the density function in the convolution model is much better than it was in the sister cell model (see Fig. 4 in [12]).

Fig. 3
Histogram (experimental data) and density function estimated by maximum-likehood technique.

7 Conclusion

We have estimated the cell cycle models that have a mixture of deterministic and probabilistic dynamics. We proposed to estimate a distribution of cell birth masses with truncated convolution of Gaussian and exponential distributions. It was shown that our convolution model fits well both the Monte Carlo simulated data and experimental data (mouse L). It is very important to be able to calculate the probability that a cell mass exceeds some level. The ability to control a pathological growth has a big significance in the biology. The saddle-point approximation method presented here is very useful in calculating such tail probabilities. This method is good not only in our convolution model, but in many other possible models, since it gives good approximations to small tail probabilities.