Modeling Microtubule Catastrophe
“The four stages of acceptance:
- This is worthless nonsense.
- This is an interesting, but perverse, point of view.
- This is true, but quite unimportant.
- I always said so.”
– JBS Haldane, 1963
The structure of the eukaryotic cell arises from the dance of the cytoskeleton, a network of protein filaments that shrink and grow dynamically. A major constituent of this network is tubulin, a monomeric protein which can polymerize to form large strands called microtubules. Along with their associated motor proteins, microtubules serve as intracellular highways for shuttling and organizing the contents of the cytoplasm. They are characterized by a slow growth phase followed by sudden rapid depolymerization known as catastrophe. A precise, quantitative understanding of microtubule dynamics would help elucidate the fundamental principles controlling growth and form inside the crowded cell.
In a 2011 publication, Gardner et al. studied the occurence of microtubule catastrophe by watching the behavior of tubulin in vitro. Here, we explore a subset of their dataset by validating the results of a control experiment, and subsequently compare two statistical models for the dynamics of microtubule catastrophe.
Fluorescently-labeled GTP-tubulin at a concentration of 12 μM was assembled onto stabilized seeds of GMPCPP-tubulin on a coverslip. Total internal reflection fluorescence (TIRF) microscopy was then used to track microtubule dynamics and record catastrophe times (Fig 2A, green). To control for phototoxicity, these measurements were repeated with unlabeled GTP-tubulin using differential interference contrast (DIC) microscopy (Fig 2A, blue). These data can be found in File 1.
Gardner et al. Fig 2A: Empirical CDFs of catastrophe times from control experiments with labeled tubulin (green, measured by TIRF microscopy) and unlabeled tubulin (blue, measured by DIC microscopy).
Using TIRF, catastrophe times were then measured for a range of labeled tubulin concentrations from 7 to 14 μM (Fig 3D). These data can be found in File 2.
Gardner et al. Fig 3D: ECDFs of catastrophe times for various concentrations of free labeled GTP-tubulin. Red: 7 μM; yellow: 9 μM; green: 10 μM; blue: 12 μM; purple: 14 μM.
All code used in the following analysis is available on the corresponding page of this website. In addition, the relevant function calls are noted where relevant, e.g. alongside the figure they were used to generate.
To begin, we set out to validate the dataset and check whether labeled tubulin catastrophizes significantly differently than unlabeled tubulin, as was done in Figure 2A. We generated ECDFs from the control data with 95% confidence intervals, confirming that they recapitulate the results of Figure 2A (plot_ECDF_controls()):
From visual inspection, it appears that the mean time to catastrophe is identically distributed for both tubulin types. We further quantified this by constructing a null hypothesis significance test (NHST) using the two-sample Kolmogorov-Smirnov test, for the null hypothesis that the distributions of mean time to catastrophe for labeled and unlabeled tubulin are identical (control_NHST()). This yielded a KS-statistic value of 0.103, and a p-value of 0.451, meaning we cannot reasonably reject the null hypothesis. Therefore, our NHST also suggests that the distributions of catastrophe times are indeed identical.
With this baseline validation complete, we proceeded to compare two potential models for microtubule catastrophe.
A Tale of Two Models
In their paper, Gardner et al. model microtubule catastrophe times with a Gamma distribution, with the following justification: we know that the Gamma distribution describes the waiting time for $\alpha$ arrivals of a Poisson process. Inside the cell, microtubule instability arises as a function of various stochastic molecular events, many of which can be thought of as Poisson processes. For example, the loss of a GTP cap from a tubulin subunit happens infrequently and at some average rate. Assuming that the events leading to catastrophe (e.g. the loss of a critical number of GTP caps) all occur at the same rate, then the time to catastrophe is the sum of exponentially distributed waiting times, which indeed makes it Gamma.1
Alternatively, we might posit that two successive biochemical processes (for example, the loss of two GTP caps) might have to occur in order to trigger catastrophe. If each process is Poisson, then this system could be modeled as the sum of two Exponential waiting times. When the rates are unique, this model is defined analytically by the following PDF:
where each $\beta$ is the arrival rate for the corresponding Poisson process. We note that when $\beta_1 = \beta_2$, this successive Poisson model collapses back into the Gamma model.
The Gamma model is parametrized by $\alpha$, the number of arrivals, and $\beta$, the arrival rate. We obtained maximum likelihood estimates (MLEs) for both parameters, along with bootstrapped 95% confidence intervals (n=10,000, parametric method, report_gamma_ci()), based on the control data for labeled tubulin:
- MLE of $\alpha$ w/ 95% CI, labeled tubulin: 2.4076, [2.0384, 2.9255]
- MLE of $\beta$ w/ 95% CI, labeled tubulin: 0.0055, [0.0045, 0.0068]
By using these parameter estimates to parametrize a Gamma distribution, we can then compare our model to the dataset (gamma_compare_MLE_ECDF()):
We observe excellent agreement between the Gamma model and the observations from experiment.
We similarly computed MLEs for the successive Poisson process model (report_diffpoisson()). CIs were not computed due to poor work ethic.
- MLE of $\beta_1$, labeled tubulin: 0.0050
- MLE of $\beta_2$, labeled tubulin: 0.0042
It is worth mentioning that we have encountered the problem of non-identifiability; it is impossible to tell from these data alone which physical processes have which of these parameters, since the results would be identical if we switched their labels. Nonetheless, we can once again parametrize our distribution using these estimates and compare it to the ECDF (diffpoisson_compare_MLE_ECDF()):
The agreement between the model and experiment is still quite good, albeit with some divergence at shorter timescales.
We can superimpose both CDFs at once to compare them more easily (compare_model_CDFs()):
It appears that the Gamma model is a better fit. We can more rigorously confirm this by calculating the Akaike information criterion (AIC) for both models; the model with the larger weight is a more likely match to the underlying generative distribution (aic()):
- Gamma weight: 0.882
- Successive Poisson weight: 0.118
The Gamma model has a much higher AIC weight than the Successive Poisson model. Therefore, it does a better job of capturing the underlying distribution of the data.
Catastrophe as a Function of Tubulin Levels
With the Gamma model emerging triumphant, we next explored the effect of titrating free tubulin concentrations on catastrophe times, using the data provided in File 2. We generated ECDFs with 95% confidence intervals, recapitulating the results of Figure 3D (plot_ECDF_exp()):
Using numerical optimization, we then obtained MLE parameter estimates for $\alpha$ and $\beta$ across all tubulin concentrations (report_gamma_params()):
|Free tubulin concentration (uM)||$\alpha$ MLE||$\beta$ MLE|
We note that $\alpha$ generally increases with concentration, but $\beta$ peaks around 10 μM before decreasing. Since we know that tubulin polymerizes more quickly at higher free tubulin concentrations, it seems that more tubulin arrival events can occur before the onset of catastrophe in this regime. However, the trend in $\beta$ may indicate that at very high concentrations, there is some stoichiometric exclusion of other necessary factors which are also needed for microtubule stability.
We showed that the in vitro dynamics of fluorescently-labeled tubulin are identical to those of unlabeled tubulin, and that a Gamma distribution can serve as an effective model for the mean time to microtubule catastrophe. Future experimental work could further dissect the role of additional molecular factors on microtubule stability, and provide additional observations with which to “assault” the Gamma model, in order to sharpen its predictive power.
We’d like to thank Professor Justin Bois and the intrepid BE/Bi 103a TAs of 2021 for an excellent and challenging quarter, as well as Gardner et al. for sharing their dataset!
Footnotes & Miscellania
It would be interesting to investigate whether the various assumptions associated with this model are valid in practice. For example, the mechanism of GTP hydrolysis may mean that different caps are lost at different rates. The “memoryless” property associated with Poisson processes may also be violated if a binding/unbinding event biases the probabilities of subsequent events, which happens in some thermodynamic models of allostery. ↩