Optimal Dynamic Instability of Microtubules

Charles S. Peskin^{1}

Abstract. Microtubules are polymers that play many important structural and functional roles within biological cells, including the separation of newly replicated chromosomes into the daughter cells during cell division. In order to catch the chromosomes that they must transport, microtubules grow out of the centrosome in each of the daughter cells. For any particular microtubule, epochs of steady growth are punctuated by episodes of rapid decay; this is known as dynamic instability. It allows for multiple attempts on the part of each microtubule to hit the small target at the center of each chromosome known as the kinetochore, where the microtubule can attach and apply traction to the chromosome. The optimal design of dynamic instability is the subject of this paper.

1991 Mathematics Subject Classification: 92C05, 92C40, 92C45

Keywords and Phrases: dynamic instability, microtubule, chromosome

1 Introduction: How to Catch and Transport a Chromosome

During cell division, the newly replicated chromosomes are pulled into the daughter cells by microtubules. This activity is organized by the centrosomes, one in each of the daughter cells, which form the two poles of the familiar mitotic spindle. Microtubules are polymers, made of protein subunits known as tubulin, that grow radially outward from the centrosomes. Dynamic instability, discovered by Mitchison and Kirschner [1], is a phenomenon concerning the assembly and disassembly of microtubules. Specifically, the individual steps of addition and removal of tubulin subunits to and from the end of a given microtubule, although random, are far from independent. Indeed, the microtubule acts like a two-state device, with a steadily growing state and a rapidly decaying state. Transitions between these states occur much more rarely than the individual steps of addition or removal of subunits.

As has been emphasized by Hill and Chen [2], dynamic instability drastically alters the statistical properties of microtubules, in comparison to the properties that would be expected on the basis of independent addition and removal of subunits. In this paper, we shall continue the exploration of this theme, from a somewhat different perspective, that of optimal design. Specifically, we shall state and solve an optimization problem that explains why dynamic instability is needed and determines certain relationships between the rate constants that characterize the assembly and disassembly of microtubules.

1Supported by the National Science Foundation (USA) under DMS/FD 92-20719. Thanks also to George Oster and David McQueen for their help in connection with this work.

Documenta Mathematica · Extra Volume ICM 1998 · III · 633–642

Mathematical theories and computer simulations of the dynamic instability of microtubules may be divided into two broad categories. First, there are the theories that simplify the microtubule by treating it as a one-dimensional polymer. Among such works are [2, 3] and also the present paper adopts this simplified point of view, which is amenable to analysis. Another possibility is to take into account the two-dimensional tubular lattice in which the subunits of the microtubule are actually arranged. This has been done in [4, 5, 6, 9]. The two-dimensional lattice models have been studied by Monte-Carlo simulation.

For a detailed review of the role of microtubules in chromosome transport, including but not limited to dynamic instability, see [7]. Dynamic instability makes possible the trial-and-error process that leads to chromosome capture by microtubules. Following capture, traction on the chromosome is generated by depolymerization of the microtubules [8].

2 Polymerization and Depolymerization of Microtubules

A typical microtubule consists of 13 protofilaments, each of which runs in a straight line, parallel to the axis of the microtubule. We shall simplify the description of the microtubule by regarding it as a one-dimensional polymer; this polymer may be thought of as representing any one of the 13 protofilaments, even though this ignores significant interactions between neighboring protofilaments, interactions which tend to coordinate their assembly and disassembly.

The subunits of a microtubule are tubulin dimers, here denoted by the symbol T. Each such tubulin dimer has two possible states, denoted T.GDP or T.GTP according to whether a guanosine diphosphate (GDP) or a guanosine triphosphate (GTP) molecule is bound to the tubulin dimer. Following the lateral cap hypothesis of Bayley et al. [5], in the simplified form appropriate to our one-dimensional model, we assume that only T.GTP can be added to a microtubule, and only T.GDP can exist in the interior of a protofilament (i.e., not at its end). Note that these rules allow the terminal subunit to be either T.GDP or T.GTP. This one bit of information will determine whether the model microtubule is in a polymerizing mode (with T.GTP at the tip), or in a depolymerizing mode (with T.GDP at the tip).

In case the terminal subunit is T.GTP, then the following polymerization reaction, driven by GTP hydrolysis (GTP → GDP + Pi, where Pi denotes inorganic phosphate) can occur

(T.GDP)_{n−1}(T.GTP) + (T.GTP) → (T.GDP)_{n}(T.GTP) + Pi  (1)

Note that the protofilament has grown by the addition of one tubulin dimer and that it still has a T.GTP subunit at its tip, so the process described by Eq.1 may be repeated indefinitely.

If, on the other hand, the terminal subunit is T.GDP, then this subunit can spontaneously dissociate:

(T.GDP)_{n} → (T.GDP)_{n−1} + (T.GDP)  (2)

and this depolymerization process, too, may be repeated indefinitely (until the microtubule has shrunk to zero length).

Conversion in either direction between the polymerizing mode (Eq.1) and the depolymerizing mode (Eq.2) may occur through the following reversible reaction, which is supposed to be rare (i.e., slow) in comparison to the reactions described by Eqs.1 and 2, above:

(T.GDP)_{n−1} + (T.GTP) ↔ (T.GDP)_{n−1}(T.GTP)  (3)

The forward reaction in Eq.3 switches the protofilament from the depolymerizing to the polymerizing mode, and the reverse reaction accomplishes the opposite. There are other possible ways to switch modes (involving phosphorylation or dephosphorylation of the terminal T.GDP or T.GTP, respectively), but we shall adhere to Eq.3 as the switching mechanism throughout this paper.

The following diagram summarizes the kinetic scheme for the assembly and disassembly of microtubules that is used in this paper:

α  α  α  α  α
A_{1}  →  A_{2}  → ⋯ →  A_{n}  → A_{n+1} → ⋯
α′ ↑↓ β′  α′ ↑↓ β′  α′ ↑↓ β′  α′ ↑↓ β′
B_{0}  ←  B_{1}  ← ⋯ ← B_{n−1} ←  B_{n}  ← ⋯
β  β  β  β  β

(4)

In this diagram, the symbol A denotes the polymerizing mode and B denotes the depolymerizing mode. The subscript on A or B denotes the total number of subunits in the polymer. Thus, A_{n} = (T.GDP)_{n−1}(T.GTP), n ≥ 1; and B_{n} = (T.GDP)_{n}, n ≥ 0. Note that B_{0} is the fixed anchor, or seed, located within the centrosome, from which the microtubule grows. Implicit in the kinetic scheme Eq.4 is the assumption that this seed has the same properties as a T.GDP molecule.

The rate constants, with dimensions of inverse time, that appear in the foregoing scheme, are defined as follows: α = rate constant for fast (Eq.1) polymerization; β = rate constant for fast (Eq.2) depolymerization; α′ = rate constant for slow (Eq.3) polymerization; β′ = rate constant for slow (Eq.3) depolymerization. Note that the polymerizing rate constants are proportional to the concentration of T.GTP in solution: α = a [T.GTP]; α′ = a′ [T.GTP], but that the depolymerizing rate constants are independent of concentration.

Implicit in our whole discussion of the microtubule as a two-state device, with a polymerizing state and a depolymerizing state, are the inequalities β′ < α and α′ < β, so that the microtubule takes many steps of polymerization with rate constant α before losing its T.GTP cap, and many steps of depolymerization with rate constant β before regaining that cap.

For the sake of comparison, however, it is also of interest to consider the special case α = α′, β = β′. This corresponds to the situation in which polymerization and depolymerization proceed without regard to the distinction between T.GTP and T.GDP, and there is no phenomenon of dynamic instability. We shall try to understand why nature does not proceed in this simple manner.

The differential equations describing an ensemble of protofilaments can now be written down by inspection of the kinetic scheme (Eq.4). Let p_{n}(t) be the probability of finding the system in state A_{n} at time t, and let q_{n}(t) be the corresponding probability for the state B_{n}. Then

\frac{dq_{0}}{dt} = β′p_{1} + βq_{1} − α′q_{0}  (5)

\frac{dp_{1}}{dt} = α′q_{0} − (α + β′)p_{1}  (6)

and for n ≥ 1

\frac{dq_{n}}{dt} = β′p_{n+1} + βq_{n+1} − (α′ + β)q_{n}  (7)
\frac{dp_{n+1}}{dt} = αp_{n} + α′q_{n} − (α + β′)p_{n+1}  (8)

Finally, the q_{n} and p_{n} are normalized according to

∑_{n=0}^{∞} q_{n} + ∑_{n=1}^{∞} p_{n} = 1  (9)

It follows from Eqs.5-8 that

∑_{k=0}^{n−1}(\frac{dq_{k}}{dt} + \frac{dp_{k+1}}{dt}) = −(αp_{n} − βq_{n})  (10)

for n ≥ 1. This will be useful in constructing a steady-state solution of Eqs.5-9.

3 Steady-State Solution [3]

In the steady state (dp_{n}/dt = dq_{n}/dt = 0 for all n), Eq.10 becomes αp_{n} = βq_{n}. Thus, we may set u_{n} = αp_{n} = βq_{n}, n ≥ 1. It then follows from the steady-state form of Eq.7 or 8 that u_{n+1} = ru_{n}, n ≥ 1, where r = (1 + α′/β)/(1 + β′/α). Thus, we have a normalizable solution if and only if r < 1, which (since all of the rate constants are positive) is equivalent to

0 < (ββ′ − αα′)  (11)

This means that depolymerization is dominant over polymerization. If the inequality Eq.11 is not satisfied, then the microtubule just grows forever and there is no steady state. From now on we shall assume that this important inequality is indeed satisfied.

According to the foregoing, the u_{n} form a geometric sequence for n ≥ 1. It is then straightforward to express all of the p_{n} and q_{n} in terms of u_{1}, and to determine u_{1} with the help of the normalization condition, Eq.9, thus completing the steady-state solution. We omit the details, but just give the following useful result:

Let N be the random variable which is the number of subunits in a protofilament, and let E[ ] denote the expected value (ensemble average) of the enclosed quantity. Then

n_{p} = E[N|N > 0] = \frac{∑_{n=1}^{∞}(p_{n} + q_{n})n}{∑_{n=1}^{∞}(p_{n} + q_{n})} = \frac{E[N]}{1 − q_{0}} = \frac{β(α + β′)}{ββ′ − αα′}  (12)

Note that n_{p} measures the average length (in subunits) of microtubules by averaging only over actual microtubules, i.e., by not including microtubules of zero length in the average.

We now compare two special cases. First suppose α = α′ and β = β′. This is the above-mentioned case in which the kinetics are indifferent to the distinction between T.GDP and T.GTP. In this case, we find n_{p} = 1/(1 − α/β) Clearly, to achieve a microtubule of any significant length (e.g., n_{p} = 100) in this situation, (α/β) must be very close to 1. On the other hand, it is also required that (α/β) < 1, or the steady-state solution does not exist. This implies that the parameters must be poised on the edge of disaster in order for the system to function!

Now consider instead the limiting case β → ∞, with α, α′, and β′ all finite. In this limit, n_{p} → (α/β′) + 1, the steady-state solution always exists, and we can make the microtubules as long as we like by choosing (α/β′) large. This is much better! The limiting case β → ∞ has other virtues as well. These will appear below.

4 Mean and Variance of the Cycle Time

In order to participate in chromosome transport, a microtubule must first grow until it hits the kinetochore of a chromosome. This being an unlikely event, repeated trials are needed. To the extent that microtubules grow in straight lines, a new trial cannot be said to begin until the microtubule shrinks all the way down to zero length and then starts to grow again, at a possibly different angle. Thus, an important random variable is the cycle time, S_{c}, which we define as the elapsed time between successive departures from the state B_{0}, in which the microtubule has zero length, see Eq.4.

The cycle time S_{c} has two components:

S_{c} = S_{0} + S_{p}  (13)

where S_{0} is the waiting time in state B_{0}, and S_{p} is the time elapsed between a given departure from B_{0} and the subsequent arrival at B_{0}. Since S_{0} and S_{p} are independent random variables, we have τ_{c} = E[S_{c}] = E[S_{0}] + E[S_{p}] and v_{c} = Var[S_{c}] = Var[S_{0}] + Var[S_{p}], where Var[ ] denotes the variance of the enclosed random variable.

On general principles concerning chemical reactions, we know that the waiting time S_{0} in the state B_{0} is exponentially distributed with mean 1/α′. It follows that E[S_{0}] = 1/α′ and Var[S_{0}] = (1/α′)^{2} Thus, to evaluate τ_{c} and v_{c}, we just need to know the mean and variance of the random variable S_{p}. These are found by assuming that the system starts (t = 0) in the state A_{1} and by treating B_{0} as an absorbing state. The reaction scheme is the same as Eq.4, except that the transition B_{0} → A_{1} is omitted:

α  α  α  α  α
A_{1} →  A_{2}  → ⋯ →  A_{n}  → A_{n+1} → ⋯
↓ β′  α′ ↑↓ β′  α′ ↑↓ β′  α′ ↑↓ β′
B_{0} ←  B_{1}  ← ⋯ ← B_{n−1} ←  B_{n}  ← ⋯
β  β  β  β  β

(14)

The differential equations are

\frac{dq_{0}}{dt} = β′p_{1} + βq_{1}  (15)

\frac{dp_{1}}{dt} = −(α + β′)p_{1}  (16)

and for n ≥ 1, we have, as before, Eqs.7 and 8. Finally, the initial conditions are p_{1}(0) = 1 with all of the other p_{n} and all of the q_{n} equal to zero at t = 0.

If we can solve this initial-value problem, then we shall have the probability density function of S_{p}, denoted ρ_{p}(t), which is given by

ρ_{p}(t) = \frac{dq_{0}}{dt} = βq_{1}(t) + β′p_{1}(t)  (17)

The initial-value problem stated above can indeed be solved in terms of Laplace transforms. Instead of inverting the Laplace transform to find ρ_{p}(t), however, we shall be content with finding the mean and variance of S_{p}, which can be evaluated directly from the Laplace transform itself. Specifically, if we define ρ̂_{p}(λ) = ∫_{0}^{∞} ρ_{p}(t) exp(−λt)dt (and similarly for all other functions of t), then we have

E[S_{p}] = ∫_{0}^{∞} tρ_{p}(t)dt = −\frac{dρ̂_{p}}{dλ} (0)  (18)

and similarly,

Var[S_{p}] = E[S_{p}^{2}] − (E[S_{p}])^{2} = \frac{d^{2}ρ̂_{p}}{dλ^{2}} (0) − (\frac{dρ̂_{p}}{dλ} (0))^{2}  (19)

In terms of the transformed variables q̂_{n} and p̂_{n}, the initial value problem becomes

λq̂_{0} − βq̂_{1} − β′p̂_{1} = 0  (20)
λp̂_{1} + (α + β′)p̂_{1} = 1  (21)

and for n ≥ 1:

(λ + α′ + β)q̂_{n} − β′p̂_{n+1} − βq̂_{n+1} = 0  (22)
(λ + α + β′)p̂_{n+1} − αp̂_{n} − α′q̂_{n} = 0  (23)

Now Eq.21 gives p̂_{1} directly, and we look for a solution of Eqs.22-23 of the following form

p̂_{n}(λ) = p̂_{1}(λ)z^{n−1}, n ≥ 1  (24)
q̂_{n}(λ) = q̂_{1}(λ)z^{n−1}, n ≥ 1  (25)

where we must require |z| < 1. With this assumed form, Eqs.22 and 23 reduce to the homogeneous 2 × 2 system

( \begin{matrix}λ_{1}z − α & −α′ \\ −β′z & λ_{2} − βz \end{matrix} ) ( \begin{matrix}p̂_{1}(λ) \\ q̂_{1}(λ) \end{matrix} ) = 0  (26)

where λ_{1} = λ + α + β′ and λ_{2} = λ + α′ + β. Of course, z is chosen so that Eq.26 has nontrivial solutions and |z| < 1. Then q̂_{1} can be found from Eq.26, since p̂_{1} is already known from Eq.21. The details are left as a (lengthy) exercise for the reader. The results, after adding the expectations of S_{p} and S_{0}, and similarly after adding their variances, are as follows:

τ_{c} = E[S_{c}] = \frac{1}{α′} + \frac{α + β}{ββ′ − αα′}  (27)

v_{c} = Var[S_{c}]
= (\frac{1}{α′} )^{2} + \frac{2β}{(α + β′)(ββ′ − αα′)} (1 + \frac{αβ′(α + β′ + α′ + β)^{2}}{(ββ′ − αα′)^{2}} )

− (\frac{α + β}{ββ′ − αα′} )^{2}  (28)

5 Optimal Design of Dynamic Instability

We are now ready to state the optimization problem that is the main subject of this paper. In order that the stochastic process of dynamic instability should proceed as regularly as possible, let us choose α, α′, β, and β′ to minimize the variance v_{c} of the cycle time, subject to given values of the mean cycle time τ_{c} and the mean length n_{p} of nonzero length microtubules.

Since there are 4 variables and 2 constraints, it should be possible to reduce the number of independent variables to 2. A convenient choice of independent variables is α′ and β. From the constraints, Eqs.12 and 27, we find

α = \frac{β(n_{p} − 1)}{Q}  (29)

β′ = \frac{α′n_{p} + β}{Q}  (30)

where

Q = (α′τ_{c} − 1)(n_{p} + \frac{β}{α′}) − α′τ_{c}(n_{p} − 1)  (31)

We can use these results to express v_{c} as a function of α′ and β only (though of course it will also contain as parameters the given constants n_{p} and τ_{c}):

v_{c}(α′, β; n_{p}, τ_{c}) = (\frac{1}{α′} )^{2} + τ_{c}^{2} (1 + 2(n_{p} − 1)\frac{α′}{β} ) − 2\frac{τ_{c}}{α′}
+ (\frac{1}{α′ + β} ) ((\frac{1}{α′} ) (2n_{p} − 1 + \frac{β}{α′} ) − 4(n_{p} − 1)τ_{c}) (32)

Our goal is to minimize v_{c} with respect to α′ and β. Let us first consider

\frac{∂v_{c}}{∂β} = −\frac{2(n_{p} − 1)}{α′(α′ + β)^{2}} ((\frac{α′ + β}{β} )^{2} (α′τ_{c})^{2} − 2(α′τ_{c}) + 1)  (33)

Since (α′ + β)/β > 1, and since n_{p} > 1, it is evident that ∂v_{c}/∂β < 0, and there can be no minimum at finite β We can, however, look for a minimum at β = ∞. Letting β → ∞, we find

v_{c}(α′, ∞; n_{p}, τ_{c}) = (\frac{1}{α′} )^{2} + (τ_{c} − \frac{1}{α′} )^{2}  (34)

which is minimized by setting

α′ = \frac{2}{τ_{c}}  (35)

To complete the solution, we need only find α and β′. Taking the limit β → ∞ in Eqs.29-31, we find

α = \frac{n_{p} − 1}{τ_{c} − \frac{1}{α′}} = \frac{2}{τ_{c}} (n_{p} − 1)  (36)

β′ = \frac{1}{τ_{c} − \frac{1}{α′}} = \frac{2}{τ_{c}}  (37)

Thus, in summary, the optimal solution is given by α′ = β′ = 2/τ_{c}; α = (n_{p} − 1)(2/τ_{c}); and β = ∞. The variance in the mean cycle time obtained in this way is given by v_{c}^{min} = τ_{c}^{2}/2, which is half that of an exponentially distributed random variable with the same mean. Note that v_{c}^{min} is independent of n_{p}.

To appreciate better the optimal solution, let us contrast it with the case obtained by setting α′ = α and β′ = β. As discussed above, this means that there is no distinction between the T.GDP and the T.GTP subunit. Under these (degenerate) circumstances, we have, after some algebra,

v_{c} = (\frac{1}{α})^{2} + \frac{2β}{(β − α)^{3}} − (\frac{1}{β − α})^{2} = (\frac{τ_{c}}{n_{p}} )^{2} (1 + (2n_{p} − 1)(n_{p} − 1)^{2}) (38)

which is asymptotic to 2n_{p}τ_{c}^{2} as n_{p} → ∞.

Thus, in the absence of a mechanism that distinguishes T.GDP from T.GTP, we find that the variance of the cycle time is a large multiple of the square of the mean cycle time, instead of being fixed at τ_{c}^{2}/2 as in the optimal solution. Such Figure 1: Near-optimal dynamic instability (above) and no dynamic instability (below); sample trajectories obtained by Monte-Carlo simulation. The horizontal axis measures time in units of τ_{c}, and the vertical axis is polymer length expressed in terms of the number of subunits. Note the extreme difference in statistical character of the trajectories, even though both have the same mean polymer length n_{p} = 25 and the same mean cycle time τ_{c} = 1. The (nearly) optimal case has many cycles of comparable duration, whereas the degenerate case has a few long cycles and a great many cycles that are much too short to be effective. This difference in statistics, which is already quite dramatic, can be tremendously accentuated by increasing the mean polymer length n_{p}. a large multiple indicates a long-tailed distribution of cycle times. Under these conditions, a microtubule that missed its target (and that would be most of them, after the first try) might spend a long time wandering up and down in length before shrinking to zero length to try again. In the case of the optimal solution, though, the cycle time is rather tightly controlled, and its variance is independent of the mean length of the microtubules. The length can therefore be made large without paying a price in terms of the variability of the cycle time. The degenerate case and a near-optimal case (finite but large β) are further contrasted in Figure 1.

References

[1] Mitchison T and Kirschner M: Dynamic instability of microtubule growth. Nature 312:237-242, 1984

[2] Hill TL and Chen Y: Phase changes at the end of a microtubule with a GTP cap. PNAS 81:5772-5776, 1984

[3] Hill TL: Introductory analysis of the GTP-cap phase change kinetics at the end of a microtubule. PNAS 81:6728-6732, 1894

[4] Chen Y and Hill TL: Monte Carlo study of the GTP cap in a five-start helix model of a microtubule. PNAS 82:1131-1135, 1985

[5] Bayley PM, Schilstra MJ, and Martin SR: Microtubule dynamic instability: numerical simulation of microtubule transition properties using a Lateral Cap model. J Cell Science 95:33-48, 1990

[6] Martin SR, Schilstra MJ, and Bayley PM: Dynamic Instability of Microtubules: Monte Carlo Simulation and Application to Different Types of Microtubule Lattice. Biophysical J 65:578-596, 1993

[7] Inoué S and Salmon ED: Force Generation by Microtubule Assembly/Disassembly in Mitosis and Related Movements. Molecular Biology of the Cell 6:1619-1640, 1995

[8] Peskin CS and Oster GF: Force production by depolymerizing microtubules: load-velocity curves and run-pause statistics. Biophysical J 69:2268-2276, 1995

[9] Tao YC and Peskin CS: Simulating the role of microtubules in depolymerization-driven transport – A Monte Carlo approach. Biophysical J. (in press)

Charles S. Peskin Courant Institute New York University 251 Mercer Street New York, NY 10012, USA peskin@cims.nyu.edu