Team:BNU-China/Model/microtubule
Microtubule Model
Introduction
To achieve our goal of arming yeasts with microtubules, we plan to display tubulin monomers on the yeast surface and culture the yeasts in a buffer that is rich in tubulin monomers under appropriate conditions to allow spontaneous microtubule polymerization on the yeast surface.
One of the key dynamics in this process is the competition between the polymerization in the solution and that on the yeast surface. We obtain the tubulins in the solution from our secretion module, and they are a valuable resource, so we want to use them as efficiently as possible. Intuitively, we know that due to the dynamical instability of microtubules, as we increase the concentration of tubulins in the solution, the number of microtubules assembled on the yeast surface will increase at first, but will gradually reach a plateau where the additional tubulins added to the system mostly polymerize in the solution. We hope to verify this intuition, find where this plateau begins and obtain the minimal concentration of tubulins in the solution that satisfies our needs.
Another thing we are interested in is how some microtubule-stabilizing agents like Taxol® can affect the polymerization process. They inhibit microtubule depolymerization by curbing GTP hydrolysis (Elie-Caille et al., 2007; Mitra and Sept, 2008). A detailed analysis of this effect can help possible future efforts to use our system to screen these agents.
In order to address the above problems, we need to be able to describe and simulate the microtubule assembly process. Taking note that stochasticity plays an important role in the polymerization and depolymerization of microtubules, we decided to construct a stochastic chemical reaction model (Mcquarrie, 1967) and conduct simulations with Monte Carlo method.
Model Formulation
We shall first describe the background settings and some basic assumptions of our model.
We consider a solution of tubulin dimers with a certain initial concentration. The tubulin dimers have already combined with GTP, which is essential for polymerization to happen. The enzymes and ions needed for the polymerization process are all optimal and held constant.
We are using a relatively complex model to study the dynamical behaviors of a microtubule population in a closed system. To simplify our calculations, we view a microtubule as a one-dimensional molecular chain, since the fact that it is composed of 13 protofilaments (Amos and Klug, 1974) does not have an essential impact on the formulation of our model. However, later we will use this information along with the fact that each tubulin dimer is about $8 \ \text{nm}$ long (Chrétien et al., 1995; Burns, 1995) to transfer the number of tubulin dimers in a microtubule to its length.
Polymerization can take place at both ends (the "+" end and the "-" end) of the microtubule. While only tubulin dimers that are combined with GTP (denoted Tu-GTP) can assemble onto the end of a microtubule, both Tu-GTP and Tu-GDP can depolymerize into the solution, although with different rate constants. The model also considers other processes involved, such as nucleation, GTP hydrolysis on microtubules and the regeneration of Tu-GTP in the solution.
Using a stochastic approach to chemical kinetics (Mcquarrie, 1967), we can capture the above dynamical processes in our model. With the parameters given by Walker et al. (1988), we can do Monte Carlo simulations to study the problems we are interested in.
Next, we shall give a detailed description of how our model is constructed. This part is folded for the convenience of those who are not interested in the technical details. Click on the arrow below to view the complete text.
Experimental data typically takes the form of rate constants, which are macroscopic descriptors of chemical processes. To do simulations on the individual microtubules, we have to transfer these data to a microscopic level. We shall now explain how this transformation can be done.
Chemical reactions first require the collision between molecules, and this collision can be viewed as a Markovian process. Therefore, chemical reactions themselves are Markovian processes. For instance, let us look at a simple reaction between two molecules:
\[ A+B{\rightarrow}C \]
From the theory of molecule collisions, we can define $\pi_{AB}$ as the probability that molecule $A$ and $B$ collide and react in a unit of time. This means that in volume $V$, the probability that $X_A$ $A$ molecules and $X_B$ $B$ molecules (that are uniformly mixed) react due to collision in time $\Delta_t$ is
\[ \pi_{AB}X_AX_B\Delta_t \]
where $\pi_{AB}$ can be seen as a microscopic rate constant of the reaction. $\pi_{AB}$ is the microscopic foundation of macroscopic deterministic rate equations, so it must be closely related to the rate constant that appears in rate equations. For the reaction described above, we have
\[ k_{AB}=\frac{\pi_{AB}V\langle X_AX_B\rangle}{\langle X_A\rangle\langle X_B\rangle} \]
Here $\langle \cdots \rangle$ stands for the ensemble average. However, in the deterministic descriptions of macroscopic dynamics, the average of a product is the same as the product of the averages, i.e. $\langle X_AX_B\rangle\approx\langle X_A\rangle\langle X_B\rangle$. Therefore, the above equation can be simplified to
\[ k_{AB}=\pi_{AB}V \]
Here, the volume $V$ appears because we use volume concentrations in the macroscopic form and number of molecules in the system in the microscopic form.
When the reaction occurs between the same kind of molecules, the distinguishable molecule pairs are
\[ \frac{X_A(X_A-1)}{2!}\approx1/2X_A*X_A \]
which yields
\[ k_{AA}\approx1/2\pi_{AA}V \]
Generalizing the above discussions, we have $k_A=\pi_A$ for first-order reactions and $k_{AB}=\pi_{AB}V, k_{AA}\approx1/2\pi_{AA}V$ for second-order reactions.
This shows that although $k$ and $\pi$ only differ by a constant, they reflect the fundamental difference between macroscopic and microscopic perspectives. It also provides the theoretical basis of comparing the results of our Monte Carlo simulation with actual experimental results.
Now that we know how to transfer between macroscopic ($k$) and microscopic ($\pi$) descriptions, we only need to analyze the microtubule dynamics in the form of $k$.
A table is provided for all the parameters in the model at the end of this section for your reference.
The dynamical model of microtubule growth can be summarized into the following steps.
\[ N=mt_t+mt_d \] \[ m_0=m_1+m_2+m_d+N \]
The meanings of the variables are given in the table below.
| Variable | Description |
|---|---|
| $N$ | The number of dimers on the microtubule |
| $mt_t$ | The number of Tu-GTP on the microtubule |
| $mt_d$ | The number of Tu-GDP on the microtubule |
| $m_0$ | The total number of tubulin dimers |
| $m_1$ | The number of Tu-GTP in the solution that have not been assembled to microtubules |
| $m_2$ | The number of Tu—GTP on the yeast surface that have not been assembled to microtubules |
| $m_d$ | The number of Tu-GDP that have depolymerized from microtubules |
The main reactions in microtubule growth are as follows. For nucleation, we have
\[ TuGTP+TuGTP\xrightarrow{k_n}(TuGTP)_2 \]
where $k_n$ is the rate constant of the nucleation reaction. The growth of a nucleus can be seen as a reversible reaction:
\[ (TuGTP)_n+TuGTP\xrightarrow{k_s}(TuGTP)_{n+1} \] \[ (TuGTP)_n+TuGTP\xleftarrow{k_p}(TuGTP)_{n+1} \]
where $k_s$ and $k_p$ are the corresponding rate constants.
A nucleus has to grow to include up to 13 tubulin dimers before it can become a relatively stable microtubule. If we denote a microtubule of length $n$ as $M_n$, we have
\[ 13 \ TuGTP\xrightarrow{}M_{13} \]
Next, we shall look at the growth of microtubules after the above process. As stated at the beginning of this section, only Tu-GTP molecules can assemble onto the end of a microtubule. For the “+” end, there are two different cases.
\[ M_nTuGTP+TuGTP\xrightarrow{k_{ptp}}M_{n+1}TuGTP \] \[ M_nTuGDP+TuGTP\xrightarrow{k_{pdp}}M_{n}TuGDP{\sim}TuGTP \]
where $k_{ptp}$ and $k_{pdp}$ are the corresponding polymerization rate constants for GTP cap and GDP cap. The “-” end is similar to the “+” end, and we denote its corresponding rate constants for GTP cap and GDP cap as $k_{stm}$ and $k_{sdm}$.
Both Tu-GTP and Tu-GDP can depolymerize from microtubules. For the “+” end we have
\[ M_nTuGTP\xrightarrow{k_{stp}}M_n+TuGTP \] \[ M_nTuGDP\xrightarrow{k_{sdp}}M_n+TuGDP \]
where $k_{stp}$ and $k_{sdp}$ are the corresponding rate constants for GTP cap and GDP cap. The “-” end is similar to the “+” end, and we denote its corresponding rate constants for GTP cap and GDP cap as $k_{stm}$ and $k_{sdm}$.
For the hydrolysis of GTP, we use an uncoupled hydrolysis model. Hydrolysis is a first-order reaction, and if we assume that only the Tu-GTP on microtubules can hydrolyze, and that all Tu-GTP on microtubules have the same hydrolysis rate constant $k_h$, we have
\[ TuGTP\xrightarrow{k_{h}}TuGDP+Pi+H^+ \]
Tu-GTP can regenerate in the solution, which is considered to be the main cause of the periodic fluctuations of the lengths of microtubules. The transformation from Tu-GDP to Tu-GTP is a relatively slow process. We describe this process with a first-order reaction with rate constant $k_{dt}$:
\[ TuGDP\xrightarrow{k_{dt}}TuGTP \]
Thus all of the major reactions in the microtubule polymerization process have been considered. Recall that the macroscopic rate constants ($k$) in the model above can be transformed into microscopic probabilities ($\pi$). With the experimental data and estimations given by Walker et al. (1998) in Table 2 below, we can run Monte Carlo simulations. The microtubules on the yeast surface are labeled differently as those in the solution, which allows us to track them separately.
| Parameter | Description | Value | Source |
|---|---|---|---|
| $k_n$ | The rate constant of microtubule nucleation | $5 \ {\mu}M^{-1}s^{-1}$ | Estimated |
| $k_p$ | The rate constant of nucleus growth | $5 \ s^{-1}$ | Estimated |
| $k_s$ | The rate constant of the reverse reaction of nucleus growth | $5 \ s^{-1}$ | Estimated |
| $k_{ptp}$ | The rate constant of polymerization at a “+” end with GTP cap | $8.9 \ {\mu}M^{-1}s^{-1}$ | Estimated |
| $k_{pdp}$ | The rate constant of polymerization at a “+” end with GDP cap | $7.5 \ {\mu}M^{-1}s^{-1}$ | Walker et al. 1998 |
| $k_{ptm}$ | The rate constant of polymerization at a “-” end with GTP cap | $4.7 \ {\mu}M^{-1}s^{-1}$ | Walker et al. 1998 |
| $k_{pdm}$ | The rate constant of polymerization at a “-” end with GDP cap | $4.0 \ {\mu}M^{-1}s^{-1}$ | Walker et al. 1998 |
| $k_{stp}$ | The rate constant of depolymerization at the “+” end with GTP cap | $44 \ s^{-1}$ | Walker et al. 1998 |
| $k_{sdp}$ | The rate constant of depolymerization at the “+” end with GDP cap | $899 \ s^{-1}$ | Walker et al. 1998 |
| $k_{stm}$ | The rate constant of depolymerization at the “-” end with GTP cap | $23 \ s^{-1}$ | Walker et al. 1998 |
| $k_{sdm}$ | The rate constant of depolymerization at the “-” end with GDP cap | $700 \ s^{-1}$ | Walker et al. 1998 |
| $k_h$ | The rate constant of GTP hydrolysis | Varies from $0$ to $5 \ s^{-1}$ in the simplified model; varies from $0$ to $10 \ s^{-1}$ in the complete model | Estimated |
| $k_{dt}$ | The rate constant of GTP regeneration | $5 \ s^{-1}$ | Estimated |
Simulation Results
Using the model constructed in the previous section, we can run Monte Carlo simulations and obtain the terminal number of the microtubules that has been formed and their length distributions under different initial conditions. By the term “terminal” we mean that overall the system has reached a state where polymerization and depolymerization reach an equilibrium. Typically, to reach the terminal state, we have to run the simulation for a time period of 30 minutes. (This is the corresponding time of real-life experiments, not the time it takes for a computer to complete the simulation.)
Even after several simplifications made in the previous section, a simulation still takes a lot of time, with one complete simulation taking up to four and a half hours. Considering the stochasticity in the model, we have to run the simulation multiple times under each initial condition and take average over the trials to obtain meaningful results, which makes it even more time-consuming to run all the simulations we need. Therefore, we decided to first run some simulations on a simplified model to give an overview of the dynamics and verify the validity of our model before using the complete model to get the key results of interest.
Preliminary Simulation Using Simplified Model
We shall now run some preliminary simulations with a simplified model to see the effects of some parameters over the microtubule dynamics and verify the validity of our model.
To simplify our model to shorten the time of a simulation, we added the following assumptions.
- Microtubules only polymerize and depolymerize on one end. Since “+” ends and “-” ends share similar characteristics and only differ by their rate constants, this is a reasonable simplification. It can cut the calculations almost by half.
- Microtubules with GTP cap and GDP cap polymerize in the same rate. This again cuts the calculations by almost a half.
- There are no regeneration of Tu-GTP in the solution. This is not too much of an oversimplification since the regeneration process is relatively slow.
Together, these assumptions reduce the time it takes to do one simulation from 4.5 hours to just about 30 minutes, which is enough for us to conduct our preliminary simulations.
We first examined the how the initial concentration of tubulin dimers $C_m$ affects the microtubule dynamics. The results are as follows.
Figure 1 The change of terminal average length (left) and terminal number (right) of microtubules as initial tubulin dimer concentration increases. $k_{ptp}=k_{pdp}=8.9{\mu}M^{-1}s^{-1}$, $k_{stp}=44s^{-1}$, $k_{sdp}=899s^{-1}$, $k_n=5{\mu}M^{-1}s^{-1}$, $k_s=5s^{-1}$.
From Figure 1 we can see that the terminal average length and the terminal number of microtubules rise with the increase of initial tubulin dimer concentration in similar trends. This is consistent with experimental results obtained by Walker et al. (1991) and Voter et al. (1991). After $C_m=10\mu M$, both of them reach a plateau. A higher concentration of tubulin dimers does not promote further polymerization, which is as we expected, since depolymerization will eventually limit the extent of microtubule assembly. The experiment done by Walker et al. (1998) showed that at $C_m=5 pM$, the polymerization and depolymerization rate are virtually the same, and microtubules can barely grow, which is in accordance with our simulation results.
We also analyzed the relationship between terminal number and the rate of GTP hydrolysis, which is controlled by $k_h$. The result is as follows.
Figure 2 The change of terminal number of microtubules as the GTP hydrolysis rate constant increases.
$k_{ptp}=k_{pdp}=8.9{\mu}M^{-1}s^{-1}$, $k_{stp}=44s^{-1}$, $k_{sdp}=899s^{-1}$, $k_n=5{\mu}M^{-1}s^{-1}$, $C_m=20{\mu}M$.
From Figure 2 we can see that the terminal number of microtubules decreases almost linearly as GTP hydrolysis accelerates. Faster GTP hydrolysis shifts the equilibrium state towards depolymerization. Here the terminal average length was omitted because the data obtained is too noisy when the terminal number is too small.
Finally, to give a more direct intuition of the microtubule growth dynamics, we chose one microtubule filament and plotted its length during the simulation period in the graph below.
Figure 3 The length of one microtubule filament in a simulation.
$k_{ptp}=k_{pdp}=8.9{\mu}M^{-1}s^{-1}$, $k_{stp}=44s^{-1}$, $k_{sdp}=899s^{-1}$, $k_n=5{\mu}M^{-1}s^{-1}$, $C_m=20{\mu}M$.
The curve in Figure 3 is similar to those obtained by Bayley et al. (1991) using Lateral Cap model and is consistent with existing experimental observations. The microtubule shifts from polymerization and depolymerization stochastically, which reflects the dynamical instability of microtubules.
This series of preliminary simulations gave us a basic understanding of the dynamics in microtubule polymerization. As explained above, each result is in accordance with the current knowledge of microtubule growth, and thus proves the validity of our model.
Key Simulation Using the Complete Model
Now that we are have a deeper understanding of our model, we can do more complex simulations with the complete model to address the key questions brought forward at the beginning, namely the effect of the initial concentration of tubulin dimers $C_m$ and the GTP hydrolysis rate constant $k_h$ (through which some microtubule-stabilizing agents take effect) on microtubule polymerization dynamics.
Here the assumptions we added in the simplified model no longer exist, and the differentiation between the tubulin dimers in the solution and those on the yeast surface is introduced. In a word, the model is made to reflect actual experimental conditions as closely as possible.
We first looked at how the initial concentration of tubulin dimers $C_m$ affect microtubule growth. The absolute value of $C_m$ is not important in this context, since different concentrations of yeasts will yield completely different outcomes with the same $C_m$. Rather, only the ratio between the tubulin dimers anchored on the yeast surface and those in the solution determine how efficient the system is. The simulation result is presented below.
Figure 4 Microtubule polymerization with different ratios between the concentration of tubulin dimers on the yeast surface and those in the solution. “Strength” is the total number of tubulin dimers that assembled on the yeast surface. $k_{ptp}=8.9{\mu}M^{-1}s^{-1}$, $k_{pdp}=7.5{\mu}M^{-1}s^{-1}$, $k_{ptm}=4.7{\mu}M^{-1}s^{-1}$, $k_{pdm}=4.0{\mu}M^{-1}s^{-1}$, $k_{stp}=44s^{-1}$, $k_{sdp}=500s^{-1}$, $k_{stm}=23s^{-1}$, $k_{sdm}=700s^{-1}$, $k_n=5{\mu}M^{-1}s^{-1} $, $C_{m1}=20{\mu}M$.
Here, “strength” is the total number of tubulin dimers assembled on the yeast surface, which reflects how efficient the tubulin dimers are being used. This quantity can be measured by fluorescence intensity, which means that our results can be tested by experiments. We can see that as the relative concentration of tubulin dimers displayed on yeasts increases, the number of dimers assembled on yeasts first increases and then stabilizes at around the ratio of 0.4. This suggests that attempting to increase the ratio to a point beyond 0.4 is unnecessary.
Next, we set the ratio to be 0.5 (to ensure the robustness of the model) and examine the polymerization on the yeast surface with different GTP hydrolysis constant $k_h$. Below is the simulation result.
Figure 5 Microtubule polymerization with different GTP hydrolysis rate constant. “Strength” is the total number of tubulin dimers that assembled on the yeast surface. Note that the vertical coordinate is logarithmic. $k_{ptp}=8.9{\mu}M^{-1}s^{-1}$, $k_{pdp}=7.5{\mu}M^{-1}s^{-1}$, $k_{ptm}=4.7{\mu}M^{-1}s^{-1}$, $k_{pdm}=4.0{\mu}M^{-1}s^{-1}$, $k_{stp}=44s^{-1}$, $k_{sdp}=500s^{-1}$, $k_{stm}=23s^{-1}$, $k_{sdm}=700s^{-1}$, $k_n=5{\mu}M^{-1}s^{-1} $, $C_{m1}=20{\mu}M$.
Note that the vertical coordinate is logarithmic. We can see that after $k_h = 2$, the strength stabilizes at a very low value, and that the change of strength around $k_h=0$ is very dramatic, reaching up to 3 orders of magnitude. Some microtubule-stabilizing agents mainly take effect by inhibiting GTP hydrolysis. Figure 5 explains how this is possible and indicates that this effect is significant on the yeast surface as well. This result confirms that our plan to use the microtubules displayed on the yeast surface to screen microtubule-stabilizing agents will be very effective.
Comparison between Simulation Results and Experimental Results
Our teammates in wet lab were able to successfully assemble microtubules on the yeast surface and used High Resolution Transmission Electron Microscopy (HRTEM) to obtain images of the microtubules (click here to see the images). Although we did not have time to conduct enough experiments to test all aspects of our model, a simulation under the conditions in this experiment can provide a simple and direct verification of the validity of our model.
According to the conditions in the experiment, we set the initial concentration of tubulin dimers to be $1.0 \ \mu \text{M}$, the volume of the system to be $100 \ \mu \text{L}$, the total number of yeasts to be $6000$, and $k_h = 0$ (since Taxol® is present in a large amount). The simulation gave us an average microtubule length of around $400 \ \text{nm}$ and about $10^3$ microtubules on each yeast. This result is consistent with our observations under HRTEM. We saw that there were a lot of microtubules formed on each yeast. Although we could not count the exact number, one can confidently say that our estimation of around $10^3$ is reasonable. The lengths of the microtubules were mostly within $200$ ~ $400 \ \text{nm}$, which is in accordance with our simulation result.
Conclusion
Using the stochastic chemical reaction model we constructed, we conducted Monte Carlo simulations to analyze microtubule polymerization dynamics. The simulation results are in accordance with existing experimental observations, which confirms the validity of the model.
We showed that the performance of the system (measured by how many tubulin dimers are polymerized onto the yeast surface) stops improving after the initial concentration ratio of dimers on the yeast surface and in the solution reaches 0.4, which helps our teammates at wet lab determine their experimental conditions.
We also illustrated that when the GTP hydrolysis rate constant $k_h$ is smaller than 1, reducing $k_h$ can boost microtubule assembly on the yeast surface dramatically. Since some microtubule-stabilizing agents such as Taxol® inhibit depolymerization by curbing GTP hydrolysis, this result proves that the agents’ stabilizing effects will still be significant in the context of our design, which confirms that our plan to use the microtubules displayed on the yeast surface to screen microtubule-stabilizing agents will be very effective.
Finally, our simulation results proved to be consistent with the experiment conducted in wet lab, which confirmed the effectiveness of our model and gave more credence to its predictions.
Amos, Linda A., and A. Klug. "Arrangement of subunits in flagellar microtubules." Journal of cell science 14.3 (1974): 523-549.
Bayley, P. M., M. J. Schilstra, and S. R. Martin. "Microtubdive dynamic instabipty: numerical simdivation of microtubdive transition properties using a Lateral Cap model. " Journal of Cell Science 95 ( Pt 1).1(1990):33.
Burns, Roy G. "α‐, β‐, and γ‐tubulins: Sequence comparisons and structural constraints." Cytoskeleton 20.3 (1991): 181-189.
Chrétien, Denis, Stephen David Fuller, and Eric Karsenti. "Structure of growing microtubule ends: two-dimensional sheets close into tubes at variable rates." The Journal of cell biology 129.5 (1995): 1311-1328.
Elie-Caille, Céline, et al. "Straight GDP-tubulin protofilaments form in the presence of taxol." Current Biology 17.20 (2007): 1765-1770.
McQuarrie, Donald A. "Stochastic approach to chemical kinetics." Journal of applied probability 4.3 (1967): 413-478.
Mitra, Arpita, and David Sept. "Taxol allosterically alters the dynamics of the tubulin dimer and increases the flexibility of microtubules." Biophysical journal 95.7 (2008): 3252-3258.
Voter, William A., E. Timothy O'Brien, and Harold P. Erickson. "Dilution‐induced disassembly of microtubules: Relation to dynamic instability and the GTP cap." Cytoskeleton 18.1 (1991): 55-62.
Walker, R. A., et al. "Dynamic instability of individual microtubules analyzed by video light microscopy: rate constants and transition frequencies." The Journal of cell biology 107.4 (1988): 1437-1448.
Walker, R. A., N. K. Pryer, and E. D. Salmon. "Dilution of individual microtubules observed in real time in vitro: evidence that cap size is small and independent of elongation rate." The Journal of cell biology 114.1 (1991): 73-81.
Copyright © 2017 BNU-China All rights reserved.
If you like this page, you can contact us: bnu_igem@163.com
