Exercise 7.2: MLE for a ligand-gated ion channel
Caltech’s Henry Lester has been studying nicotinic acetylcholine receptors for decades. A schematic of the receptor, which is an ion channel, and its five (four unique) subunits, adapted from the book Physical Biology of the Cell 2, is shown below.
When acetylcholine ligands bind to the α units (the binding sites are shown as divots in the schematic), the channel opens and ions can flow through.
In a 1995 paper (Labarca, et al.), Labarca and coworkers in the Lester lab measured the (normalized) voltage across the ion channel for varying concentrations of acetylcholine. They did this for five different genotypes, shown below, where an asterisk denotes a mutation.
\(\alpha_2\beta\gamma\delta\) (wild type)
\(\alpha_2^*\beta\gamma\delta\)
\(\alpha_2\beta\gamma^*\delta\)
\(\alpha_2^*\beta\gamma^*\delta^*\)
\(\alpha_2\beta^*\gamma^*\delta^*\)
You can download the data set here.
We can work out theoretically that the fraction of time that an ion channel is open follows a logistic equation,
\begin{align} f_\mathrm{open}(c) = \frac{1}{1 + \mathrm{e}^{-\beta F(c)}}, \end{align}
where \(c\) is the concentration of acetylcholine, \(\beta = 1/k_B T\) is the inverse thermal energy, and \(F(c)\) is the Bohr parameter (named after Christian Bohr, the father of Nobel laureate Niels Bohr and of Danish national team footballer and mathematician Harald Bohr, and grandfather of Nobel laureate Aage Bohr), given by
\begin{align} F(c) = 2\Delta E_\mathrm{u} + 2 k_B T \ln\left(\frac{1+c/K_\mathrm{do}}{1+c/K_{dc}}\right). \end{align}
Here, \(\Delta E_\mathrm{u} = E_\mathrm{unbound}^\mathrm{closed} - E_\mathrm{unbound}^\mathrm{open}\) is the energy difference between the closed and open state of the ion channel without any acetylcholine bound, \(K_\mathrm{do}\) is the dissociation constant for acetylcholine binding the receptor when it is open, and \(K_\mathrm{dc}\) is the dissociation constant for acetylcholine binding the receptor when it is closed.
If we assume that in the limit of large concentrations the ion channel is completely open, then the normalized voltage \(v(c)\) across the channel is given by \(f_\mathrm{open}(c)\); \(v(c) = f_\mathrm{open}(c)\).
Before proceeding to the specific questions, I want to note why we are doing this problem.
It serves to review techniques of parameter estimation by optimization (in this case maximum likelihood estimation). We will use very similar techniques in a Bayesian context in the next couple of weeks.
It provides and example and a watch-out for how fragile statistical modeling is. You will need to think about the structure of your model and investigate its failure modes. You can do this with techniques you learned last term, and we will spend a lot of time with Bayesian approaches to this.
a) A common approach to parameter estimation is to find the so-called “best-fit” parameters and report them. “Best fit” is interpreted in a least squares sense and finding those parameters involves a least squares calculation in which the sum of the squares of the residuals is minimized. This approach is equivalent to performing a maximum likelihood estimate of a model in which the data are Normally distributed about the theoretical curve with a variance \(\sigma^2\) that is the same for all concentrations. Specifically,
\begin{align} v_i \mid c_i, \beta\Delta E_\mathrm{u}, K_\mathrm{do}, K_\mathrm{dc}, \sim \text{Norm}(v(c_i;\beta\Delta E_\mathrm{u}, K_\mathrm{do}, K_\mathrm{dc}), \sigma). \end{align}
For each of the five genotypes get the “best fit” parameters \(\beta\Delta E_\mathrm{u}\), \(K_\mathrm{do}\), \(K_\mathrm{dc}\). You do not need to compute confidence intervals for these estimates in this part, though you may in the next parts of the problem. You may with to refer to this lesson from last term for reference.
b) For each “best fit” parameter, plot the “best fit curve” along with the data. This is the theoretical curve parametrized by the parameters you found in part (a).
c) In my experience, the vast majority of analyses stop after part (b). TO GREAT PERIL! This model/data combination is an example of a case where quantitative interpretation of the estimates you get for the parameters would lead to erroneous conclusions. I will tell you now that this model is not identifiable with these measurements. (Remember, when we say a model is not identifiable, it means that we cannot unambiguously determine the parameter values. That is, two or more parameter sets are observationally equivalent.) You job in this part of the problem is to demonstrate that the model is not identifiable.