Building a model¶

The statistical model¶

Not all spindles will be measured to be exactly $\varphi$ µm in length. Rather, there may be some variation about $\varphi$ due to natural variation and measurement error. So, we would expect measured length of spindle i to be

$$\begin{array}{r}{l}_i=\varphi +e_i,\end{array}$$

where $e_i$ is the noise component of the $i$the datum.

So, we have a theoretical model for spindle length, $l=\varphi$, and to get a fully generative model, we need to model the errors $e_i$. A reasonable model assumes

1. Each measured spindle’s length is independent of all others.

2. The variability in measured spindle length is Normally distributed.

With these assumptions, we can write the probability density function for $l_i$ as

$$\begin{array}{r}f(l_i;\varphi ,\sigma )=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp [-\frac{(l_i-\varphi {)}^2}{2{\sigma }^2}].\end{array}$$

Since each measurement is independent, we can write the joint probability density function of the entire data set, which we will define as $\mathbf{l}=\{l_1,l_2,\dots \}$, consisting of $n$ total measurements.

$$\begin{array}{r}f(\mathbf{l};\varphi ,\sigma )=\prod _{i}f(l_i;\varphi ,\sigma )=\frac{1}{(2\pi {\sigma }^2{)}^{n/2}}\exp [-\frac{1}{2{\sigma }^2}\sum _i(l_i-\varphi {)}^2].\end{array}$$

We can write this more succinctly, and perhaps more intuitively, as

$$\begin{array}{r}{l}_i\sim \text{Norm}(\varphi ,\sigma )\mathrm{\forall }i.\end{array}$$

We will generally write our models in this format, which is easier to parse and understand. Note that in writing this generative model, we have necessarily introduced another parameter, $\sigma$, the standard deviation parametrizing the Normal distribution. So, we have two parameters in our model, $\varphi$ and $\sigma$.

The cartoon model¶

The three key principles of this “cartoon” model are:

1. The total amount of tubulin in the droplet or cell is conserved.

2. The total length of polymerized microtubules is a function of the total tubulin concentration after assembly of the spindle. This results from the balances of microtubule polymerization rate with catastrophe frequencies.

3. The density of tubulin in the spindle is independent of droplet or cell volume.

The mathematical model¶

From these principles, we need to derive a mathematical model that will provide us with testable predictions. The derivation follows below (following the derivation presented in the paper), and you may read it if you are interested. Since our main focus here is building a statistical model, you can skip ahead to to the final equation, where we define a mathematical expression relating the spindle length, $l$ to the droplet diameter, $d$, which depends on two parameters, $\gamma$ and $\varphi$. Nonetheless, it is important to see how a models such as this one is derived.

Principle 1 above (conservation of tubulin) implies

$$\begin{array}{r}{T}_0{V}_0=T_1(V_0-V_{\mathrm{s}})+T_{\mathrm{s}}{V}_{\mathrm{s}},\end{array}$$

where $V_0$ is the volume of the droplet or cell, $V_{\mathrm{s}}$ is the volume of the spindle, $T_0$ is the total tubulin concentration (polymerized or not), $T_1$ is the tubulin concentration in the cytoplasm after the the spindle has formed, and $T_{\mathrm{s}}$ is the concentration of tubulin in the spindle. If we assume the spindle does not take up much of the total volume of the droplet or cell ($V_0\gg V_{\mathrm{s}}$, which is the case as we will see when we look at the data), we have

$$\begin{array}{r}{T}_1\approx T_0-\frac{V_{\mathrm{s}}}{V_0}{T}_{\mathrm{s}}.\end{array}$$

The amount of tubulin in the spindle can we written in terms of the total length of polymerized microtubules, $L_{\mathrm{MT}}$ as

$$\begin{array}{r}{T}_s{V}_{\mathrm{s}}=\alpha L_{\mathrm{MT}},\end{array}$$

where $\alpha$ is the tubulin concentration per unit microtubule length. (We will see that it is unimportant, but from the known geometry of microtubules, $\alpha \approx 2.7$ nmol/µm.)

We now formalize assumption 2 into a mathematical expression. Microtubule length should grow with increasing $T_1$. There should also be a minimal threshold $T_{\min }$ where polymerization stops. We therefore approximate the total microtubule length as a linear function,

$$\begin{array}{r}{L}_{\mathrm{MT}}\approx \{\begin{array}{ccl}0& & T_1\le T_{\min }\\ \beta (T_1-T_{\min })& & T_1>T_{\min }.\end{array}\end{array}$$

Because spindles form in Xenopus extract, $T_0>T_{\min }$, so there exists a $T_1$ with $T_{\min }<T_1<T_0$. Thus, going forward, we are assured that $T_1>T_{\min }$. So, we have

$$\begin{array}{r}{V}_{\mathrm{s}}\approx \alpha \beta \frac{T_1-T_{\min }}{T_{\mathrm{s}}}.\end{array}$$

With insertion of our expression for $T_1$, this becomes

$$\begin{array}{r}{V}_{\mathrm{s}}\approx \alpha \beta (\frac{T_0-T_{\min }}{T_{\mathrm{s}}}-\frac{V_{\mathrm{s}}}{V_0}).\end{array}$$

Solving for $V_{\mathrm{s}}$, we have

$$\begin{array}{r}{V}_{\mathrm{s}}\approx \frac{\alpha \beta }{1+\alpha \beta /V_0}\frac{T_0-T_{\min }}{T_{\mathrm{s}}}=\frac{V_0}{1+V_0/\alpha \beta }\frac{T_0-T_{\min }}{T_{\mathrm{s}}}.\end{array}$$

We approximate the shape of the spindle as a prolate spheroid with major axis length $l$ and minor axis length $w$, giving

$$\begin{array}{r}{V}_{\mathrm{s}}=\frac{\pi }{6}l{w}^2=\frac{\pi }{6}{k}^2{l}^3,\end{array}$$

where $k\equiv w/l$ is the aspect ratio of the spindle. We can now write an expression for the spindle length as

$$\begin{array}{r}l\approx {\left(\frac{6}{\pi k^2}\frac{T_0-T_{\min }}{T_{\mathrm{s}}}\frac{V_0}{1+V_0/\alpha \beta }\right)}^{\frac{1}{3}}.\end{array}$$

For small droplets, with $V_0\ll \alpha \beta$, this becomes

$$\begin{array}{r}l\approx {\left(\frac{6}{\pi k^2}\frac{T_0-T_{\min }}{T_{\mathrm{s}}}{V}_0\right)}^{\frac{1}{3}}={\left(\frac{T_0-T_{\min }}{k^2{T}_{\mathrm{s}}}\right)}^{\frac{1}{3}}d,\end{array}$$

where $d$ is the diameter of the spherical droplet or cell. So, we expect the spindle size to increase linearly with the droplet diameter for small droplets.

For large $V_0$, the spindle size becomes independent of droplet size;

$$\begin{array}{r}l\approx {\left(\frac{6\alpha \beta }{\pi k^2}\frac{T_0-T_{\min }}{T_{\mathrm{s}}}\right)}^{\frac{1}{3}}.\end{array}$$

Indentifiability of parameters¶

We measure the microtubule length $l$ and droplet diameter $d$ directly from the images. We can also measure the spindle aspect ratio $k$ directly from the images. Thus, we have four unknown parameters, since we already know α ≈ 2.7 nmol/µm. The unknown parameters are:

We would like to determine all of these parameters. We could measure them all either in this experiment or in other experiments. We could measure the total tubulin concentration $T_0$ by doing spectroscopic or other quantitative methods on the Xenopus extract. We can $T_{\min }$ and $T_s$ might be assessed by other in vitro assays, though these parameters may by strongly dependent on the conditions of the extract.

Importantly, though, the parameters only appear in combinations with each other in our theoretical model. Specifically, we can define two parameters,

$$\begin{array}{rl}\gamma & ={\left(\frac{T_0-T_{\min }}{k^2{T}_{\mathrm{s}}}\right)}^{\frac{1}{3}}\\ \varphi & =\gamma {\left(\frac{6\alpha \beta }{\pi }\right)}^{\frac{1}{3}}.\end{array}$$

We can then rewrite the general model expression in terms of these parameters as

$$\begin{array}{r}l(d)\approx \frac{\gamma d}{{(1+(\gamma d/\varphi {)}^3)}^{\frac{1}{3}}}.\end{array}$$

If we tried to determine all four parameters from this experiment only, we would be in trouble. This experiment alone cannot distinguish all of the parameters. Rather, we can only distinguish two combinations of them, which we have defined as $\gamma$ and $\varphi$. This is an issue of identifiability. We may not be able to distinguish all parameters in a given model, and it is important to think carefully before the analysis about which ones we can identify. Even with our work so far, we are not quite done with characterizing identifiability.

Visualizing the mathematical model¶

Let’s take a quick look at the mathematical model so we can see how the curve looks. It’s best to nondimensionalize the diameter by $\varphi$, giving

$$\begin{array}{r}\frac{l(d)}{\varphi }\approx \frac{\gamma d/\varphi }{{(1+(\gamma d/\varphi {)}^3)}^{\frac{1}{3}}}.\end{array}$$

So, we will plot $l(/\varphi$ versus $d/\varphi$, which means we choose units of length to be $\varphi$.

[3]:

hv.extension("bokeh")

def theor_spindle_length(gamma, d):
    """Compute spindle length using mathematical model"""
    return gamma * d / np.cbrt(1 + (gamma * d)**3)

d = np.linspace(0, 20, 200)

def plot_theor(gamma):
    return hv.Curve(
        data=(d, theor_spindle_length(gamma, d)),
        kdims=['d/φ'],
        vdims=['l/φ'],
        label=f'γ = {gamma}',
    ).opts(
        color=hv.Cycle(list(bokeh.palettes.Blues7[1:-1][::-1]))
    )

plots = [plot_theor(gamma) for gamma in [0.03, 0.1, 0.3, 0.7, 1.0]]

hv.Overlay(plots)

[3]:

The curve grows from zero to a plateau at $l=\varphi$, more rapidly for larger $\gamma$. We can more carefully characterize the limiting behavior.

Limiting behavior¶

For large droplets, with $d\gg \varphi /\gamma$, the spindle size becomes independent of $d$, with

$$\begin{array}{r}l\approx \varphi .\end{array}$$

Conversely, for $d\ll \varphi /\gamma$, the spindle length varies approximately linearly with diameter.

$$\begin{array}{r}l(d)\approx \gamma d.\end{array}$$

Note that the expression for the linear regime gives bounds for $\gamma$. Obviously, $\gamma >0$, lest we get negative spindle lengths. Because $l\le d$, lest the spindle not fit in the droplet, we also have $\gamma \le 1$.

Importantly, if the experiment is done in the regime where $d$ is large (and we do not really know a priori how large that is since we do not know the parameters $\varphi$ and $\gamma$), we cannot tell the difference between the two models, since they are equivalent in that regime. Further, if the experiment is in this regime the model is unidentifiable because we cannot resolve $\gamma$.

This sounds kind of dire, but this is actually a convenient fact. The second model is more complex, but it has the simpler model, model 1, as a limit. Thus, the two models are in fact commensurate with each other. Knowledge of how these limits work also enhances the experimental design. We should strive for small droplets. And perhaps most importantly, if we didn’t consider the second model, we might automatically assume that droplet size has nothing to do with spindle length if we simply did the experiment in larger droplets.

Generative model¶

We have a theoretical model relating the droplet diameter to the spindle length. Let us now build a generative model. For spindle, droplet pair i, we assume

$$\begin{array}{r}{l}_i=\frac{\gamma d_i}{{(1+(\gamma d/\varphi {)}^3)}^{\frac{1}{3}}}+e_i.\end{array}$$

We will assume that $e_i$ is Normally distributed with variance ${\sigma }^2$. This leads us to our statistical model.

$$\begin{array}{rl}& {\mu }_i=\frac{\gamma d_i}{{(1+(\gamma d_i/\varphi {)}^3)}^{\frac{1}{3}}},\\ & l_i\sim \text{Norm}({\mu }_i,\sigma )\mathrm{\forall }i,\end{array}$$

which is equivalently stated as

$$\begin{array}{r}{l}_i\sim \text{Norm}(\frac{\gamma d_i}{{(1+(\gamma d_i/\varphi {)}^3)}^{\frac{1}{3}}},\sigma )\mathrm{\forall }i.\end{array}$$

Importantly, note that this model builds upon our first model. Generally, when doing modeling, it is a good idea to build more complex models on your initial baseline model such that the models are related to each other by limiting behavior. This gives you a continuum of model and a sound basis for making comparisons among models.

Note that we are assuming the droplet diameters are known. When we generate data sets for prior predictive checks, we will randomly generate them from about 20 µm to 200 µm, since this is the range achievable with the microfluidic device.