# Our numerical workhorses
import numpy as np
# Image processing tools
import skimage.io
With your new image processing skills, you can start to work on more complicated images. For this problem, you will analyze a FRAP experiment and then perform parameter estimation to determine a diffusion coefficient and a binding rate for two molecules. We will be taking a simplified approach, but there is more sophisticated analysis we can do to get better estimates for the phenomenological coefficients.
The data set comes from Nate Goehring. The images are taken of a C. elegans one-cell embryo expressing a GFP fusion to the PH domain of Protein Lipase C delta 1 (PH-PLCd1). This domain binds PIP2, a lipid enriched in the plasma membrane. By using FRAP, we can investigate the dynamics of diffusion of the PH-PLCd1/PIP2 complex on the cell membrane, as well as the binding/unbinding dynamics of PH-PLCd1.
So, if $c$ is the concentration of the PH-PLDd1/PIP2 complex on the membrane and $c_\mathrm{cyto}$ is the concentration of PH-PLCd1 in the cytoplasm (assumed to be spatially uniform since diffusion in the cytoplasm is very fast), the dynamics are described by a reaction-diffusion equation.
\begin{align} \frac{\partial c}{\partial t} = D\left(\frac{\partial^2 c}{\partial x^2} + \frac{\partial^2 c}{\partial y^2}\right) + k_\mathrm{on} c_\mathrm{cyto} - k_\mathrm{off} c. \end{align}Here, $k_\mathrm{on}$ and $k_\mathrm{off}$ are the phenomenological rate constants for binding and unbinding to PIP2 on the membrane, and $D$ is the diffusion coefficient for the PH-PLCd1/PIP2 complex on the membrane.
In their paper, the authors discuss techniques for analyzing the data taking into account the fluorescence recovery of the bleached region in time and space. For simplicity here, we will only consider recovery of the normalized mean fluorescence. If $I(t)$ is the mean fluorescence of the bleached region and $I_0$ is the mean fluorescence of the bleached region immediately before photobleaching, we have, as derived in the paper,
\begin{align} I_\mathrm{norm}(t) \equiv \frac{I(t)}{I_0} &= 1 - f_b\,\frac{4 \mathrm{e}^{-k_\mathrm{off}t}}{d_x d_y}\,\psi_x(t)\,\psi_y(t),\\[1mm] \text{where } \psi_i(t) &= \frac{d_i}{2}\,\mathrm{erf}\left(\frac{d_i}{\sqrt{4Dt}}\right) -\sqrt{\frac{D t}{\pi}}\left(1 - \mathrm{e}^{-d_i^2/4Dt}\right), \end{align}where $d_x$ and $d_y$ are the extent of the photobleached box in the $x$- and $y$-directions, $f_b$ is the fraction of fluorophores that were bleached, and $\mathrm{erf}(x)$ is the error function. Note that this function is defined such that the photobleaching event occurs at time $t = 0$.
We measure $I(t)$, $I_0$, $d_x$, and $d_y$. We can also measure $f_b$ as
\begin{align} f_b \approx 1 - \frac{I(0^+)}{I_0}, \end{align}though we will consider this a parameter to estimate. In practice, the normalized fluorescent recovery does not go all the way to unity. This is because the FRAP area is a significant portion of the membrane, and we have depleted fluorescent molecules. We should thus adjust our equation to be
\begin{align} I_\mathrm{norm}(t) \equiv \frac{I(t)}{I_0} &= f_f\left(1 - f_b\,\frac{4 \mathrm{e}^{-k_\mathrm{off}t}}{d_x d_y}\,\psi_x(t)\,\psi_y(t)\right), \end{align}where $f_f$ is the fraction of fluorescent species left. So, we have four parameters to use in regression, the physical parameters of interest, $D$ and $k_\mathrm{off}$, and $f_f$ and $f_b$.
The FRAP images come in a TIFF stack, which is a single TIFF file containing multiple frames. You can load these with the skimage.io.ImageCollection class.
# Load in TIFF stack
fname = '../data/t6_data/goehring_frap_data/PH_138_A.tif'
ic = skimage.io.ImageCollection(fname, conserve_memory=False)
# How long is it?
print('There are {0:d} frames.'.format(len(ic)))
A problem with this particular TIFF stack is that the first frame is loaded in with a NumPy array that has all frames.
print('Shape of first "frame":', ic[0].shape)
print('Do the indices match?', np.allclose(ic[0][1,:,:], ic[1]))
You can either use ic[0][0,:,:] in place of ic[0] in your code, or you can just use the the three-dimensional NumPy array.
a) Extract the mean normalized fluorescence versus time from each of the TIFF stacks for the experimental repeats. Note that important information is contained in the associated README file. You can download the data set here.
b) Perform regressions to find $D$ and $k_\mathrm{off}$ for each. Report your final estimates for $D$ and $k_\mathrm{off}$.
c) (10 pts extra credit) Repeat part (b) with a hierarchical model.