Homework 9.4: Determining dissociation constants (100 pts)

Data set 1 download, Data set 2 download

---

In this problem, we will work with data from a paper by Rasson and coworkers. The authors were investigating the biochemistry of Spire-actin interactions. Spire is an actin binding protein that can nucleate actin filaments. In particular, it has four domains (called \(S_A\), \(S_B\), \(S_C\), and \(S_D\)), which bind monomeric actin. These four domains, acting in concert, can line up actin monomers to help in nucleation. In this problem, we will determine the dissociation constant, \(K_d\), describing binding of \(S_D\) to monomeric actin.

The strategy to determine \(K_d\) is to perform a titration experiment and then use nonlinear regression to determine \(K_d\). Consider the chemical reaction describing \(S_D\)-actin binding.

\begin{align} \text{actin}\cdot S_D \rightleftharpoons \text{actin} + S_D, \end{align}

which has dissociation constant \(K_d\). (In this system actin does not polymerize and stay monometic because it is treated with Latrunculin.) Let \(c_a\) be the equilibrium concentration of actin and \(c_d\) be the equilibrium concentration of \(S_D\), and \(c_{ad}\) be the equilibrium concentration of bound actin-\(S_D\). Then, at equilibrium,

\begin{align} K_d = \frac{c_a c_d}{c_{ad}}. \end{align}

If we started with a total actin concentration of \(c_a^0\) and a total \(S_D\) concentration of \(c_d^0\), we also have

\begin{align} c_a^0 = c_a + c_{ad}, \\[1mm] c_d^0 = c_d + c_{ad}, \end{align}

by conservation of mass. With these relations, we can now write \(c_{ad}\) in terms of \(c_a^0\) and \(c_d^0\), which are known quantities (this is what we pipetted into our solution).

\begin{align} K_d &= \frac{(c_a^0 - c_{ad})(c_d^0 - c_{ad})}{c_{ad}},\\[1mm] \Rightarrow\;&\;\;c_{ad}^2 - (K_d + c_a^0 + c_d^0)c_{ad} + c_a^0 c_d^0 = 0. \end{align}

The solution to this quadratic equation gives \(c_{ad}\) as a function of \(K_d\). Note that we must choose one of the two roots, the one that is physical. The physical root satisfies \(0 < c_{ad} < \min(c_a^0, c_d^0)\). In this case, it is

\begin{align} c_{ad} = \frac{1}{2}\left(K_d + c_a^0 + c_d^0 - \sqrt{\left(K_d + c_a^0 + c_d^0\right)^2 - 4c_a^0c_d^0}\right). \end{align}

In a titration experiment, we fix \(c_d^0\) and vary \(c_a^0\), and measure \(c_{ad}\) to get a curve. From the curve, we can perform a regression to get \(K_d\).

The problem with this approach is that we do not have a direct way of measuring \(c_{ad}\). The authors instead employed fluorescence anisotropy. I will not go into the details here of how it works, but will simply say that larger complexes rotate more slowly, and therefore give a higher fluorescence anisotropy signal (which is dimensionless) than do smaller complexes.

So, the authors fluorescently tagged \(S_D\). We will call this molecule \(S_{D^*}\), with concentration \(c_{d^*}\). When free in solution, this molecule gives an anisotropy signal of \(r_f\). When bound to actin, it gives an anisotropy signal of \(r_b\). So, the total anisotropy signal we detect is

\begin{align} r = \frac{1}{c_{d^*}^0}\,\left(r_f c_{d^*} + r_b c_{ad^*}\right). \end{align}

Clearly, when all \(S_{D^*}\) is free, the anisotropy signal is \(r_f\) and when all \(S_{D^*}\) is bound to actin, the signal is \(r_b\). Remembering conservation of mass, \(c_{d^*} = c_{d^*}^0 - c_{ad^*}\), we have

\begin{align} r = \frac{1}{c_{d^*}^0}\,\left(r_f (c_{d^*}^0 - c_{ad^*}) + r_b c_{ad^*}\right) = r_f + \frac{r_b-r_f}{c_{d^*}^0}\, c_{ad^*}. \end{align}

Now, returning to our equilibrium expression, we have

\begin{align} c_{ad^*} = \frac{1}{2}\left(K_d^* + c_a^0 + c_{d^*}^0 - \sqrt{\left(K_d^* + c_a^0 + c_{d^*}^0\right)^2 - 4c_a^0c_{d^*}^0}\right), \end{align}

so we can write the measured anisotropy \(r\) as a function of \(K_d^*\) and the known quantities \(c_a^0\) and \(c_{d^*}^0\). Note that we now have three parameters in our mathematical model, \(K_d^*\), \(r_f\), and \(r_b\), since the latter two are not known a priori.

This is all fine and good, but if we do this regression, we are measuring the dissociation constant of \(S_{D^*}\), not \(S_D\). To get \(K_d\), we can use the fact that we know \(K_d^*\) from dissociation experiments described above. Now, say we add monomeric actin, \(S_{D^*}\), and \(S_D\) to a reaction mixture. Then, we have two reactions going on.

\begin{align} \text{actin-}S_D &\rightleftharpoons \text{actin} + S_D \\[1mm] \text{actin-}S_{D^*} &\rightleftharpoons \text{actin} + S_{D^*}, \end{align}

with equilibrium constants \(K_d\) and \(K_d^*\), respectively. In this case, we have five equations describing equilibrium, the two equilibrium expressions and three conservation of mass expressions.

\begin{align} K_d &= \frac{c_a c_d}{c_{ad}} \\[1mm] K_d^* &= \frac{c_a c_{d^*}}{c_{ad^*}} \\[1mm] c_a^0 &= c_a + c_{ad} + c_{ad^*}\\[1mm] c_d^0 &= c_d + c_{ad} \\[1mm] c_{d^*}^0 &= c_{d^*} + c_{ad^*}. \end{align}

These five equations can be rearranged to give

\begin{align} c_a^3 + \beta c_a^2 + \gamma c_a + \delta = 0, \end{align}

with

\begin{align} \beta &= K_d + K_d^* + c_d^0 + c_{d^*}^0 - c_a^0, \\[1mm] \gamma &= K_d(c_{d^*}^0 - c_a^0) + K_d^*(c_d^0 - c_a^0) + K_d K_d^* \\[1mm] \delta &= -K_d K_d^* c_a^0. \end{align}

So, we can solve this third order polynomial for \(c_a\). We can then compute \(c_{d^*}\) and \(c_{ad^*}\) using the equilibrium and mass conservation relations for \(S_{D^*}\) as

\begin{align} c_{d^*} &= \frac{K_d^* c_{d^*}^0}{K_d^* + c_a} \\[1mm] c_{ad^*} &= \frac{c_a c_{d^*}^0}{K_d^* + c_a}. \end{align}

Given these expressions for \(c_{ad^*}\) and \(c_{d^*}\), we can compute the anistropy as a function of \(K_d\), \(K_d^*\), and the known quantities \(c_a^0\), \(c_d^0\), and \(c_{d^*}^0\).

This looks like a complicated function for the anisotropy. This is why researchers have consistently fit competition anisotropy data with approximate (wrong) functions. In fact, the way most people have done this makes approximations that neglect the most dynamic part of the curve! In practice, though, this is not too much of a challenge; we need only to solve a cubic polynomial.

Your task in this problem is to use all of the available data to get a good estimate for \(K_d\) (and you will get estimates for other parameters as you do so). You should use a hierarchical model.

Big hint: To compute the anisotropy in the competition assay, you will need to solve the cubic equation above. To do so, you can include the following functions in the functions block of your Stan code to do so.

real solve_ca(
  real Kd,
  real Kd_star,
  real ca0,
  real cd0,
  real cd_star0
) {
  // Coefficients of cubic
  real a = 1.0;
  real b = Kd + Kd_star + cd0 + cd_star0 - ca0;
  real c = Kd * (cd_star0 - ca0) + Kd_star * (cd0 - ca0) + Kd * Kd_star;
  real d = -Kd * Kd_star * ca0;

  real f = (3.0 * c / a - b^2 / a^2) / 3.0;
  real g = (2.0 * b^3 / a^3 - 9.0 * b * c / a^2 + 27.0 * d / a) / 27.0;
  real h = g^2 / 4.0 + f^3 / 27.0;

  if (f == 0 && g == 0 && h == 0) {
    // All three roots are real and equal
    return cbrt(-d / a);
  }

  if (h <= 0) {
    // All three roots are real
    real i = sqrt(g^2 / 4.0 - h);
    real j = cbrt(i);
    real k = acos(-g / 2.0 / i);
    real m = cos(k / 3.0);
    real n = sqrt(3.0) * sin(k / 3.0);
    real p = -b / 3.0 / a;

    real x1 = 2 * j * cos(k / 3.0) - b / 3.0 / a;
    real x2 = p - j * (m + n);
    real x3 = -b / 3.0 / a;

    // Return the physical root
    if (x1 <= ca0 && x1 >= 0.0) return x1;
    else if (x2 <= ca0 && x2 >= 0.0) return x2;
    else return x3;
  }

  // If we got here, we have one real root and two complex roots
  real r = -g / 2.0 + sqrt(h);
  real s;
  if (r >= 0) s = cbrt(r);
  else s = -cbrt(-r);

  real t = -g / 2.0 - sqrt(h);
  real u;
  if (t >= 0) {
    u = cbrt(t);
  }
  else {
    u = -cbrt(-t);
  }

  // Return the real root
  return s + u - b / 3.0 / a;
}


real comp_anisotropy(
  real Kd,
  real Kd_star,
  real ca0,
  real cd0,
  real cd_star0,
  real rf,
  real rb_minus_rf
) {
  // Concentration of free actin
  real ca = solve_ca(Kd, Kd_star, ca0, cd0, cd_star0);

  // Get concentration of bound labeled SpirD
  real cad_star = ca * cd_star0 / (Kd_star + ca);

  return rf + rb_minus_rf * cad_star / cd_star0;
}