Homework 7.2: Models for microtubule catastrophe (70 pts)

Data set download

---

We have thoroughly investigated the process by which microtubules undergo catastrophe using data from the Gardner, Zanic, et al. paper. We used an Exponential, Gamma, Weibull, and our custom two-step distribution to model the catastrophe times. We have consistently found that the Gamma model works best. It does have the shortcoming,though, that it does not directly match a story that might arise from chemical kinetics, which we strongly suspect would regulate microtubule castastrophe. We expect an integer number $m$ of Poisson processes to arrive sequentially in order for catastrophe to occur. The Exponential and two-step model are special cases we have already considered where $m=1$ and $m=2$, respectively. In this problem, we will consider a model for arbitrary $m$ and assess which $m$ gives the most plausible generative model. We will again use the experiment run with a tubulin concentration of 12 µM as our data set.

The probability density function for the time $t$ to catastrophe triggered by the arrival of $m$ Poisson processes may be shown to be

$$\begin{array}{r}f(t\mid \mathit{\beta },m)=\sum _{j=1}^m\frac{{\mathrm{e}}^{-{\beta }_{j}t}}{{\beta }_j^{m-2}\prod _{k=1,k\ne j}^m\frac{{\beta }_k-{\beta }_j}{{\beta }_j{\beta }_k}}.\end{array}$$

The expression is a bit cleaner when written in terms of ${\tau }_i=1/{\beta }_i$.

$$\begin{array}{r}f(t\mid \mathit{\tau },m)=\sum _{j=1}^m\frac{{\tau }_j^{m-2}{\mathrm{e}}^{-t/{\tau }_j}}{\prod _{k=1,k\ne j}^m({\tau }_j-{\tau }_k)}.\end{array}$$

For clarity, the likelihoods for the first few $m$ are

$$\begin{array}{rl}f(t\mid {\tau }_1,1)& =\frac{{\mathrm{e}}^{-t/{\tau }_1}}{{\tau }_1},\\ f(t\mid {\tau }_1,{\tau }_2,2)& =\frac{{\mathrm{e}}^{-t/{\tau }_1}}{{\tau }_1-{\tau }_2}+\frac{{\mathrm{e}}^{-t/{\tau }_2}}{{\tau }_2-{\tau }_1}=\frac{{\mathrm{e}}^{-t/{\tau }_2}-{\mathrm{e}}^{-t/{\tau }_1}}{{\tau }_2-{\tau }_1},\\ f(t\mid {\tau }_1,{\tau }_2,{\tau }_3,3)& =\frac{{\tau }_1^{}{\mathrm{e}}^{-t/{\tau }_1}}{({\tau }_1-{\tau }_2)({\tau }_1-{\tau }_3)}+\frac{{\tau }_2{\mathrm{e}}^{-t/{\tau }_2}}{({\tau }_2-{\tau }_1)({\tau }_2-{\tau }_3)}+\frac{{\tau }_3{\mathrm{e}}^{-t/{\tau }_3}}{({\tau }_3-{\tau }_1)({\tau }_3-{\tau }_2)},\\ f(t\mid {\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4,4)& =\frac{{\tau }_1^2{\mathrm{e}}^{-t/{\tau }_1}}{({\tau }_1-{\tau }_2)({\tau }_1-{\tau }_3)({\tau }_1-{\tau }_4)}+\frac{{\tau }_2^2{\mathrm{e}}^{-t/{\tau }_2}}{({\tau }_2-{\tau }_1)({\tau }_2-{\tau }_3)({\tau }_2-{\tau }_4)}\\ & +\frac{{\tau }_3^2{\mathrm{e}}^{-t/{\tau }_3}}{({\tau }_3-{\tau }_1)({\tau }_3-{\tau }_2)({\tau }_3-{\tau }_4)}+\frac{{\tau }_4^2{\mathrm{e}}^{-t/{\tau }_4}}{({\tau }_4-{\tau }_1)({\tau }_4-{\tau }_2)({\tau }_4-{\tau }_3)}.\end{array}$$

Note that these probability distributions assume that none of the ${\tau }_j$’s are equal, and you should explicitly ensure this in your calculations (Hint: You may want to read the solutions to Homework 5.2, in particular the Stan model for the two-step model, to see how to implement this.)

a) Build a model for arbitrary $m$ and code it up in Stan. Draw samples out of this model for values of $m$ ranging from 1 to whatever you think is reasonable.

b) Compare the models for various $m$. Which value of $m$ gives the best predictive model based on a model comparison problem?

c) This is another example where I think the model comparison is unnecessary and should actually be avoided. Without directly doing a model comparison by computing a LOO or WAIC like you did in part (b), interpret the results of your sampling to advocate for a physical picture of how catastrophe proceeds.