Analysis of models of flagellar growth¶
(c) 2019 Justin Bois. With the exception of pasted graphics, where the source is noted, this work is licensed under a Creative Commons Attribution License CC-BY 4.0. All code contained herein is licensed under an MIT license.
This document was prepared at Caltech with financial support from the Donna and Benjamin M. Rosen Bioengineering Center.
This document was generated from a Jupyter notebook, which may be downloaded here.
# Our numerical workhorses
import numpy as np
import scipy.integrate
import scipy.optimize
# Import plotting tools
import bokeh.plotting
import bokeh.io
bokeh.io.output_notebook()
The balance point model¶
The balance point model gives the rate of growth of a flagellum (in dimensionless form) as
\begin{align} \frac{\mathrm{d}\ell}{\mathrm{d}t} = \frac{1}{\ell} - 1. \end{align}
Numerical integration of the balance point model¶
We can solve this numerically using the scipy.integrate.odeint() function. It uses the Hindmarsh algorithm, intelligently dealing with potential stiffness in the equations.
You can look at the documentation for scipy.integrate.odeint() here. The typical function call to scipy.integrate.odeint() is of the form
scipy.integrate.odeint(f, y0, t, args=()).
scipy.integrate.odeint() solves the system of ODEs
\begin{align} \frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t} =f(\mathbf{y}, t), \end{align}
where $\mathbf{y}$ is a vector and $f(\mathbf{y},t)$ is vector-valued. Thus, scipy.integrate.odeint() takes as its first argument a function that returns an array containing the right hand side of the system of ODEs you are computing. When you define this function, it must be of the form
f(y, t, *args)
where the *args indicates the parameters on which the function depends. The second argument to scipy.integrate.odeint() is the initial condition, again stored as an array. The third argument is an array of time points for which you want the solution to the ODEs. Finally, args is a tuple containing the other parameters to be passed into the function f.
In the present case, we have only a single ODE. We just have to code up the right hand side and solve!
def balance_point_rhs(ell, t):
"""
Right hand side of balance equation. I use `ell` for length because
`l` looks like a one.
"""
return 1 / ell - 1
Now that we have defined that function, we specify the initial condition. In our case, we will watch how the flagellum grows from a small nub of length 0.1. (Remember, $\ell = 0$ is unphysical under this model.)
ell0 = np.array([0.1])
Finally, we specify the dimensionless time points we want for our solution.
t = np.linspace(0, 4, 200)
Now, we call scipy.integrate.odeint() to do the solution.
ell = scipy.integrate.odeint(balance_point_rhs, ell0, t).flatten()
Done! That's it. Now, we just have to plot the solution.
# Set up the figure
p = bokeh.plotting.figure(plot_width=650,
plot_height=350,
y_range=[0,1],
x_axis_label='dimensionless time',
y_axis_label='dimensionless flag. length')
# Specify the glyphs
p.line(t, ell, line_width=3)
# Render
bokeh.io.show(p);
Shrinkage dynamics¶
Interestingly, we can also watch a flagellum shrink. Let's start it at length $\tilde{\ell} = 2$.
# Initial condition for shrinkage
ell0 = np.array([2])
# Solve ODE
t = np.linspace(0, 5, 200)
ell_shrink = scipy.integrate.odeint(balance_point_rhs, ell0, t).flatten()
# Plot both together
p = bokeh.plotting.figure(plot_width=650,
plot_height=350,
y_range=[0,2],
x_axis_label='dimensionless time',
y_axis_label='dimensionless flag. length')
# Populate glyphs
p.line(t, ell, line_width=3)
p.line(t, ell_shrink, line_width=3, color='orange')
bokeh.io.show(p);
The model predicts that a flagellum will reach its steady state more rapidly by growth than by shrinkage.
Comparing the balance point model to experiment¶
Wallace Marshall's lab has investigated flagellar growth extensively. In Engel, Ludington, and Marshall, J. Cell Biol., 187, 81-89, 2009, they measured the growth of a flagellum after a pH shock (which eliminates the flagella). They did the repeatedly (2026 times) and averaged the length of the flagella as a function of time after pH shock. I digitized mean data from their paper.
t_exp = np.array([ 2.48447205, 7.20496894, 12.42236025, 17.14285714,
22.11180124, 27.08074534, 31.80124224, 36.77018634,
41.73913043, 46.45962733, 51.92546584, 56.89440994,
61.86335404, 66.58385093, 71.80124224, 76.52173913,
81.49068323, 86.45962733, 91.18012422, 101.11801242,
111.05590062, 120.49689441, 130.93167702, 140.62111801,
150.31055901, 160.49689441])
flag_len_exp = np.array([ 0.12704174, 0.12704174, 0.12704174, 2.03266788,
3.81125227, 5.56442831, 6.5553539 , 7.59709619,
8.53720508, 9.19782214, 9.40108893, 9.8076225 ,
10.08711434, 10.03629764, 10.36660617, 10.44283122,
10.59528131, 10.59528131, 10.82395644, 10.95099819,
11.15426497, 11.30671506, 11.45916515, 11.38294011,
11.63702359, 11.4845735 ])
# Make plot of data and curve fit
p = bokeh.plotting.figure(plot_width=650,
plot_height=350,
x_axis_label='time (min)',
y_axis_label='flag. length (µm)')
p.circle(t_exp, flag_len_exp, size=7)
bokeh.io.show(p);
There is a lag after the pH shock, presumably for the cell to recover, and then growth begins. We will fit these data with the balance point model. We will take time zero to be the fourth data point, the first where we see a significant flagellum. The fit parameters are $\alpha$, $\beta$, and also $\ell_0$, the initial length of the flagellum when the growth kicks in (since we do not know this exactly). To perform the curve fit, we need to first write a function that will return the theoretical flagellar length as a function of time for the given parameters.
def flag_len_theor(t, alpha, beta, ell0):
"""
Compute theoretical flagellar length given parameters
"""
# Dimensionless ell0
ell0_dimless = ell0 * alpha / beta
# Dimenesionless time
t_dimless = t * alpha**2 / beta
# Solve dimensionless equations
ell = scipy.integrate.odeint(balance_point_rhs, ell0_dimless, t_dimless)
# Return with dimensional units
return ell.flatten() * beta / alpha
Now that we have this function, we just provide an initial guess for the parameter values and use scipy.optimize.curve_fit() to get the optimal parameters. We will guess, based on eyeballing the curve, that $\alpha = 0.5$ µm/s, $\beta = 4.5$ µm$^2$/s, and $\ell_0 = 2$ µm.
# Perform the curve fit with initial guess p0
p0 = (0.5, 4.5, 2)
popt, _ = scipy.optimize.curve_fit(flag_len_theor,
t_exp[3:] - t_exp[3],
flag_len_exp[3:],
p0=p0)
# Print best-fit parameter values
print("""
Best fit parameter values
α: {0:.2f} µm/min
β: {1:.2f} µm^2/min
l_0: {2:.2f} µm
""".format(*popt))
We can qualitatively look at how well the model fits the data by plotting.
# Compute theoretical curve
t = np.linspace(0, 160, 200)
flag_len_model = flag_len_theor(t, *popt)
# Add curve fit to the plot
p.line(t + t_exp[3], flag_len_model, line_width=3, color='dimgray')
bokeh.io.show(p);
By eye, this is convincing support for the balance point model. However, as described in the lecture notes, this model does not capture the interdependencies of the growth of the two flagella.
Caveat (a major caveat)¶
The procedure we have just done amounts to a maximumum likelihood estimate of the parameters of the theoretical model based on the data. Implicit in this estimate is the generative statistical model. It assumed that every data point is Normally distributed about the theoretical curve, and further that these Normal distributions are independent and have the same variance. There is a whole world of hurt you could be opening yourself up to in taking this naive approach. I encourage you to take my course on data analysis to learn more about better ways to perform this kind of analysis.
The adjusted balance point model¶
The adjusted balance point model takes into account cytoplasmic concentrations of precursor for flagellum building. The dynamical equations are
\begin{align} \frac{\mathrm{d}\ell_1}{\mathrm{d}t} &= \beta\, \frac{n(t-\ell_1/v_a)}{\ell_1} - \alpha,\\[1em] \frac{\mathrm{d}\ell_2}{\mathrm{d}t} &= \beta \, \frac{n(t-\ell_2/v_a)}{\ell_2} - \alpha,\\[1em] \frac{\mathrm{d}n}{\mathrm{d}t} &= -\beta n\left(\frac{1}{\ell_1} + \frac{1}{\ell_2}\right) + 2\alpha. \end{align}
Here, $v_a$ is the speed of anterograde motion of IFT trains and $n$ is the number of precursor units available in the nearby cytoplasm. We nondimensionalize using the definitions
\begin{align} \gamma &= \frac{\beta}{\alpha}, \;\; u = \frac{v_a}{\alpha}, \;\; \ell_1 = \gamma n_\mathrm{tot}\,\tilde{\ell_1}, \\[1em] \ell_2 &= \gamma n_\mathrm{tot}\,\tilde{\ell_2},\;\; t = \frac{\beta n_\mathrm{tot}}{\alpha^2}\,\tilde{t},\;\; n = n_\mathrm{tot}\tilde{n}, \end{align}
giving dimensionless dynamical equations
\begin{align} \frac{\mathrm{d}\ell_1}{\mathrm{d}t} &= \frac{n(t-\ell_1/u)}{\ell_1} - 1,\\[1em] \frac{\mathrm{d}\ell_2}{\mathrm{d}t} &= \frac{n(t-\ell_2/u)}{\ell_2} - 1,\\[1em] \gamma^{-1}\,\frac{\mathrm{d}n}{\mathrm{d}t} &= - n(t)\left(\frac{1}{\ell_1} + \frac{1}{\ell_2}\right) + 2, \end{align}
where we have dropped the tildes for notational convenience. This is the set of delay differential equations we have to solve. We will solve them using Euler stepping. We first write a function for the time derivatives of $\ell_1$, $\ell_2$, and $n$.
def dx_dt(x, n_1, n_2, gamma, u):
"""
Right hand side of DDE for shared pool balance point model.
x = ell_1, ell_2, n
n_1 = n(t - ell_1/u)
n_2 = n(t - ell_2/u)
"""
return np.array([n_1 / x[0] - 1,
n_2 / x[1] - 1,
gamma * (2 - x[2] * (1/x[0] + 1/x[1]))])
Next, we write a function to solve the adjusted balance point DDE using Euler stepping. We have to specify some points in the "past" because of the delay/.
def solve_balance_point_dde(x0, gamma, u, dt=0.001, t_stop=5):
"""
Use Euler integration to solve ODEs
"""
# Make sure initial lengths are physical
if gamma * (x0[0] + x0[1]) + x0[2] > 1:
raise RuntimeError('gamma * (ell_1 + ell_2) + n must be <= 1')
# Starting time point (in the past)
t0 = -max(x0[0] / u, x0[1] / u)
# Time points to consider
t_points = np.concatenate(((t0,), np.arange(0, t_stop, dt)))
# Initialize output array
x = np.empty((len(t_points), 3))
x[0:2,:] = np.array([x0, x0])
# Do Euler stepping
for i, t in enumerate(t_points[2:]):
# Do linear interpolation of past points to find past n
n_1 = np.interp(t - x[i+1,0]/u, t_points[:i+2], x[:i+2,2])
n_2 = np.interp(t - x[i+1,1]/u, t_points[:i+2], x[:i+2,2])
# Take Euler step
x[i+2] = x[i+1] + dt * dx_dt(x[i+1], n_1, n_2, gamma, u)
return t_points, x
Now, let's see if the adjusted balance point model can fit growth of a single microtubule (assuming, by symmetry, that the two microtubules grow identially in the pH shock experiment). We will take $v_a = 120$ µm/min, as measured by Kaminsky, et al., PNAS, 90, 5519–5523, 1993.
First, we write a function to generate theoretical values of the growth of the pair of microtubules in the pH shock experiment for specified time points.
v_a = 120 # µm/min
def flag_len_adj_theor(t, alpha, beta, n_tot, ell_0):
"""
Theoretical flagellar length.
"""
# Parameters
gamma = beta / alpha
u = v_a / alpha
# Dimensionless ell_0
ell_0_dimless = ell_0 / gamma / n_tot
# Dimensionless time
t_dimless = t * alpha**2 / beta / n_tot
# Solve dimensionless equations
t_solve, x = solve_balance_point_dde(
np.array([ell_0_dimless, ell_0_dimless, 1 - 2*gamma*ell_0_dimless]),
gamma,
u,
dt=0.0001,
t_stop=t_dimless.max())
# Interpolate to values we want
ell_dimless = np.interp(t_dimless, t_solve, x[:,0])
# Convert to real lengths
return ell_dimless * gamma * n_tot
Now that we have a theoretical function, we can perform the curve fit.
# Curve fit
p0 = (0.1, 0.1, 30, 2)
popt, _ = scipy.optimize.curve_fit(flag_len_adj_theor, t_exp[3:] - t_exp[3],
flag_len_exp[3:], p0=p0)
# Compute useful quantities from results
alpha, beta, n_tot, _ = popt
gamma = beta / alpha
u = v_a / alpha
# Steady states
n_ss = alpha * n_tot / (alpha + 2 * beta)
ell_ss = gamma * n_ss
# Show results
print("""
Best fit parameter values
α: {0:.3f} µm/min
β: {1:.3f} µm/min
n_tot: {2:.2f} µm
l_0: {3:.2f} µm
""".format(*popt))
Let's plot the result and see how we did.
# Compute theoretical curve
t = np.linspace(0, 160, 200)
flag_len_model = flag_len_adj_theor(t, *popt)
# Make plot of solution
p.line(t + t_exp[3], flag_len_model, line_width=3, color='tomato')
bokeh.io.show(p);
Not surprisingly, the more complicated model also fits the data well. In fact, the models are almost exactly coincident, as seem by the overlap of the two theoretical lines.
Now, let's see how the system will respond to a severing experiment. We will take as an initial condition one flagellum being at the steady state length and the other being close to zero. The total amount of tubulin is as in before, minus what was in the severed flagellum.
new_n_tot = n_tot - ell_ss
t, x = solve_balance_point_dde(np.array([ell_ss / new_n_tot, 0.01, n_ss/new_n_tot]),
gamma,
u,
dt=0.001,
t_stop=1)
Now that we have the solution, we can plot the result. The severed filament is in red and the unsevered filament is in blue.
# Make plot of solution
p = bokeh.plotting.figure(plot_width=650,
plot_height=350,
x_axis_label='time (min)',
y_axis_label='flag. length (µm)')
t_plot = t * beta / alpha**2 * new_n_tot
p.line(t_plot, x[:,0] * gamma * new_n_tot, line_width=3)
p.line(t_plot, x[:,1] * gamma * new_n_tot, line_width=3, color='crimson')
bokeh.io.show(p);