Competition anisotropy

In this case study, we demonstrate how the low-level interface of EQTK can be integrated into other code that is just-in-time compiled using Numba. Specifically, we will perform a Markov chain Monte Carlo (MCMC) calculation of competition anisotropy titration curves to infer dissociation constants.

As our example system, we will consider determination of a dissociation constant between two biomolecules using a competition fluorescence anisotropy. The technique relies on the fact that a fluorophore that is excited with polarized light undergoes rotational diffusion during the time it is fluorescing, which leads to depolarized emitted light. The fluorescence anisotropy is related to the rotational diffusion coefficient of the fluorophore. An unbound fluorophore has faster rotational diffusion, and therefore lower anisotropy, than a fluorophore that is bound to another molecule.

If F is the fluorescent molecule and A is its binding partner (our molecule of interest), we have the chemical reaction

\begin{align} \mathrm{AF} \rightleftharpoons \mathrm{A} + \mathrm{F}. \end{align}

The anisotropy, \(r\), is given by

\begin{align} r = \frac{1}{c_\mathrm{F}^0}(r_\mathrm{f} c_{\mathrm{F}} + r_\mathrm{b} c_{\mathrm{AF}}), \end{align}

the concentration-weighted average of the anisotropy of free fluorophore, \(r_\mathrm{f}\), and that of bound fluorophore, \(r_\mathrm{b}\). We define \(K_\mathrm{d,f}\) as the dissociation constant for this reaction. By adjusting the concentration of A and measuring the anisotropy, we can determine \(K_\mathrm{d,f}\) by performing a regression.

With the dissociation constant of the fluorophore in hand, we can then consider the coupled chemical reactions

\begin{align} &\mathrm{AF} \rightleftharpoons \mathrm{A} + \mathrm{F},\\[1em] &\mathrm{AB} \rightleftharpoons \mathrm{A} + \mathrm{B}, \end{align}

where the binding of A and B is of interest. This assay is called a competition assay, as the molecule of interest, B, competes with the fluorophore, F, for binding of A. We wish to determine the dissociation constant, \(K_\mathrm{d}\), of the second reaction. To do this, the concentration of A and F are fixed and the concentration of B is varied. We then perform a regression to get \(K_\mathrm{d}\). This involves a coupled equilibrium, and EQTK can solve for \(c_\mathrm{F}\) and \(c_\mathrm{AF}\), which are required to compute the theoretical anisotropy, \(r\).

Before we begin our analysis, let’s import the necessary packages, including HoloViews for plotting, also setting some nice defaults for the plots.

import numpy as np
import pandas as pd

import eqtk

import numba

import holoviews as hv

hv.extension("bokeh")
hv.opts.defaults(
    hv.opts.Scatter(
        alpha=0.75,
        color=hv.Cycle(["#4c78a8", "#f58518", "#e45756"]),
        fontsize=dict(labels=12, legend=8, title=12),
        height=250,
        padding=0.05,
        show_grid=True,
        size=5,
        toolbar="above",
        width=350,
    )
)

WARNING:param.main: pandas could not register all extension types imports failed with the following error: cannot import name 'ABCIndexClass' from 'pandas.core.dtypes.generic' (/Users/bois/opt/anaconda3/lib/python3.8/site-packages/pandas/core/dtypes/generic.py)

The system to analyze and data

We will analyze data from Rasson, et al. (J. Mol. Biol., 2015). The aim is to determine the dissociation coefficient between a domain of an actin nucleator, termed SD, and actin monomers. Let’s take a look at the anisotropy data from the two experiments, one to measure fluorophore-actin binding, and the other a competition assay to measure SD-actin binding.

First let’s load in the data and take a look at the data format. (If you want to run this notebook, the data set is available here).

df = pd.read_csv('anisotropy_data.csv', comment='#')

df.head()

	[SD] (µM)	[F] (µM)	[actin] (µM)	anisotropy	trial
0	0.0	0.005	31.10	115.2	dissociation
1	0.0	0.005	16.55	114.0	dissociation
2	0.0	0.005	8.28	113.9	dissociation
3	0.0	0.005	4.14	112.7	dissociation
4	0.0	0.005	2.07	109.0	dissociation

The initial concentrations of SD, the fluorophore F, and actin, as well as the measured anisotropy are given for dissociation and competition assays. Let’s plot these data.

diss_plot = hv.Scatter(
    data=df.loc[df['trial']=='dissociation', :],
    kdims='[actin] (µM)',
    vdims='anisotropy',
    label='fluorophore dissociation',
)

comp_plot = hv.Scatter(
    data=df.loc[df['trial']=='competition', :],
    kdims='[SD] (µM)',
    vdims='anisotropy',
    label='competition',
)

diss_plot + comp_plot

Statistical model

For our analysis using Markov chain Monte Carlo, we will assume that the experimental anisotropy curves follow the theoretical curves above with Normally distributed homoscedastic error. The value of \(K_\mathrm{d,f}\) does not vary between the dissociation and competition experiments, but the measured anisotropy of free and bound fluorophore can vary. As a result, we have seven parameters to determine.

Parameter	Description
\(K_\mathrm{d}\)	dissociation constant for SD and actin
\(K_\mathrm{d,f}\)	dissociation constant for fluorophore and actin
\(r_\mathrm{f,d}\)	anisotropy of free F in dissociation experiment
\(r_\mathrm{b,d}\)	anisotropy of bound F in dissociation experiment
\(r_\mathrm{f,c}\)	anisotropy of free F in competition experiment
\(r_\mathrm{b,c}\)	anisotropy of bound F in competition experiment
\(\sigma\)	homoscedastic error in measurement

We define by \(\mathbf{c}_i^0\) the initial concentrations of all species for anisotropy measurement \(i\). The ordering of the species is

\begin{align} \mathbf{c}_i^0 = (c_\mathrm{A, i}^0, c_\mathrm{B,i}^0, c_\mathrm{AB,i}^0, c_\mathrm{AF,i}^0)^\mathsf{T}. \end{align}

We define \(r_i\) to be the anisotropy found in measurement \(i\). Finally, we define by \(\hat{r}_i\) to be the theoretical anisotropy as computed by solving the coupled equilibrium problem and then computing

\begin{align} \hat{r}_i = \frac{1}{c_\mathrm{F,i}^0}(r_f c_{\mathrm{F,i}} + r_b c_{\mathrm{AF,i}}). \end{align}

We will choose Log-Normal priors for the dissociation constants (equivalent to choosing Normal priors for the free energies of the species), Normal priors for the bound and free anisotropy parameters, and a Half-Normal prior for the error \(\sigma\). Assuming dimensionless dissociation constants, the statistical model is then

\begin{align} &\ln K_\mathrm{d} \sim \text{Norm}(0, 10), \\[1em] &\ln K_\mathrm{d,f} \sim \text{Norm}(0, 10), \\[1em] &r_\mathrm{f,d} \sim \text{Norm}(100, 30), \\[1em] &r_\mathrm{b,d} \sim \text{Norm}(100, 30), \\[1em] &r_\mathrm{f,c} \sim \text{Norm}(100, 30), \\[1em] &r_\mathrm{b,c} \sim \text{Norm}(100, 30), \\[1em] &\sigma \sim \text{HalfNorm}(0, 10),\\[1em] &r_i \sim \text{Norm}(\hat{r}_i, \sigma). \end{align}

MCMC sampling

To perform the sampling, we will use a Metropolis algorithm. (Note: This is not a preferred sampling algorithm; Hamiltonian Monte Carlo, for example, would perform better. We use Metropolis here to demonstrate how EQTK’s capabilities can be integrated into a larger calculation.) We will not show the details of the sampling algorithm here, but the source code is shown at the end of this document. To be concise, will run the .py file containing the sampling functions so we have access to them here.

%run metropolis_hastings.py

Let’s look at the doc string of mh_sample(), the main sampling function.

help(mh_sample)

Help on function mh_sample in module __main__:

mh_sample(logpost, x0, sigma=None, discrete=None, args=(), n_burn=1000, n_steps=1000, tune_interval=100, variable_names=None, return_acceptance_rate=False)
    Parameters
    ----------
    logpost : function
        The function to compute the log posterior. It has call
        signature `logpost(x, *args)`, and must by compiled with
        numba in nopython mode.
    x0 : ndarray, shape (n_variables,)
        The starting location of a walker in parameter space.
    sigma : ndarray, shape (n_variables,)
        The standard deviations for the proposal distribution.
        If None, takes all values to be one.
    discrete : ndarray of bools, shape (n_variables,)
        discrete[i] is True if variable i is discrete and False
        otherwise. If None (default), all entries are False.
    args : tuple
        Additional arguments passed to `logpost()` function.
    n_burn : int, default 1000
        Number of burn-in steps.
    n_steps : int, default 1000
        Number of steps to take after burn-in.
    tune_interval : int, default 100
        Number of steps to use when determining acceptance
        fraction for tuning.
    variable_names : list, length n_variables
        List of names of variables. If None, then variable names
        are sequential integers.
    return_acceptance_rate : bool, default False
        If True, also return acceptance rate.

    Returns
    -------
    output : DataFrame
        The first `n_variables` columns contain the samples.
        Additionally, column 'lnprob' has the log posterior value
        at each sample.

Importantly, the sampling function, mh_sample() requires an inputted log posterior function that is JITted with Numba. This means that we should use EQTK’s low-level interface to compute the equilibrium concentrations and ultimately the anisotropy.

Numba’d anisotropy

We will use eqtk.solveNK(), a JITted low-level function for computing the equilibrium concentrations. EQTK’s low-level functions require that the input be as contiguous Numpy arrays in C order. Further, instead of specifying the equilibrium constants, \(\mathbf{K}\), we supply their natural logarithm. Finally, the low-level functions return the natural logarithm of the dimensionless concentrations. With these considerations in mind, and also remembering the ordering of the concentrations (actin, actin-SD, F, actin-SD, actin-F), we can write a JITted function to compute the anisotropy.

@numba.njit
def anisotropy(x0, N, logKd, logKdf, rf, rb):
    """User EQTK to compute the anisotropy."""
    logx = eqtk.solveNK(x0, N, np.array([logKd, logKdf]))
    x = np.exp(logx)

    return (rf * x[:,2] + rb * x[:,4]) / x0[:,2]

Specifying the log posterior

We can now specify the log posterior according to our statistical model. Note that we call the anisotropy() function in the log likelihood function. The anisotropy() function in turn calls eqtk.solveNK(). Because all of these functions are JITted with Numba, the entire log posterior can also be JITted.

@numba.njit
def log_prior(logKd, logKdf, rf_d, rf_c, rb_d, rb_c, sigma):
    """Log prior, enforcing positivity constraints."""
    if rf_d <= 0 or rf_c <= 0 or rb_d <= 0 or rb_c <= 0 or sigma <= 0:
        return -np.inf

    lp = -(logKd ** 2) / 200
    lp -= logKdf ** 2 / 200
    lp -= (rf_d - 100) ** 2 / 1800
    lp -= (rf_c - 100) ** 2 / 1800
    lp -= (rb_d - 100) ** 2 / 1800
    lp -= (rb_c - 100) ** 2 / 1800
    lp -= sigma ** 2 / 200

    return lp


@numba.njit
def log_likelihood(
    logKd, logKdf, rf_d, rf_c, rb_d, rb_c, sigma, x0_d, x0_c, r_d, r_c, N
):
    """Log likelihood, with anisotropies computed using EQTK."""
    r_d_theor = anisotropy(x0_d, N, logKd, logKdf, rf_d, rb_d)
    r_c_theor = anisotropy(x0_c, N, logKd, logKdf, rf_c, rb_c)

    log_like = -np.sum((r_d - r_d_theor) ** 2) / (2 * sigma ** 2)
    log_like -= np.sum((r_c - r_c_theor) ** 2) / (2 * sigma ** 2)
    log_like -= (len(r_d) + len(r_c)) * np.log(sigma)

    return log_like


@numba.njit
def log_posterior(params, x0_d, x0_c, r_d, r_c, N):
    """Log posterior with ordering of parameters:
    logKd, logKdf, rf_d, rf_c, rb_d, rb_c, sigma
    """
    logKd, logKdf, rf_d, rf_c, rb_d, rb_c, sigma = params

    lp = log_prior(logKd, logKdf, rf_d, rf_c, rb_d, rb_c, sigma)

    if lp == -np.inf:
        return -np.inf

    log_like = log_likelihood(
        logKd, logKdf, rf_d, rf_c, rb_d, rb_c, sigma, x0_d, x0_c, r_d, r_c, N
    )

    return lp + log_like

Setting up inputs

We now need to set up the inputs for the sampler. Everything must be Numpy arrays, so we will start by extracting the measured anisotropies and initial concentrations from the data frame containing the experimental measurements.

# Convert data to Numpy array
r_d = df.loc[df["trial"] == "dissociation", "anisotropy"].values
r_c = df.loc[df["trial"] == "competition", "anisotropy"].values
c0_d = df.loc[
    df["trial"] == "dissociation", ["[actin] (µM)", "[SD] (µM)", "[F] (µM)"]
].values
c0_c = df.loc[
    df["trial"] == "competition", ["[actin] (µM)", "[SD] (µM)", "[F] (µM)"]
].values

EQTK requires initial concentrations of all species, so we need to add in zero concentrations for the bound species, actin-SD and actin-F.

# Add in the zero initial concentrations of bound species
c0_d = np.concatenate((c0_d, np.zeros((len(c0_d), 2))), axis=1)
c0_c = np.concatenate((c0_c, np.zeros((len(c0_c), 2))), axis=1)

Finally, the low-level functions work in dimensionless units, so we need to nondimensionalize the concentrations.

# Make dimensionless
solvent_density = eqtk.water_density(T=293.15, units='µM')
x0_d = c0_d / solvent_density
x0_c = c0_c / solvent_density

The equilibrium constants are explored by the sampler, so we do not need to specify those, but we do need to make the stoichiometric matrix. Note that it must be of float data type to be used in the low-level EQTK functions.

# Stoichiometric matrix
N = np.array([[1, 1, 0, -1, 0], [1, 0, 1, 0, -1]], dtype=float)

Sampling

We are almost ready to perform the sampling. We first need to tell the sampler where to start.

params_0 = np.array(
    [
        0,   # log Kd
        0,   # log Kdf
        100, # rf_d
        100, # rf_c
        60,  # rb_d
        60,  # rb_c
        10   # sigma
    ],
    dtype=float)

For convenience of output, we should also specify the names of the variables.

variable_names = ("logKd", "logKdf", "rf_d", "rf_c", "rb_d", "rb_c", "sigma")

Now we can sample! After sampling, we can compute the dissociation constants by exponentiating their logarithms and converting to units of micromolar.

df_samples = mh_sample(
    log_posterior,
    params_0,
    args=(x0_d, x0_c, r_d, r_c, N),
    n_burn=10000,
    n_steps=10000,
    variable_names=variable_names,
)

# Convert to dimensional Kd's
df_samples["Kd"] = np.exp(df_samples["logKd"]) * solvent_density
df_samples["Kdf"] = np.exp(df_samples["logKdf"]) * solvent_density

Visualizing the result

We can now visualize the result, plotting the posterior distribution of the two parameters of interest (the dissociation constants).

points = hv.Points(
    data=df_samples.iloc[::10, :],
    kdims=[('Kd', 'Kd (µM)'), ('Kdf', 'Kd* (µM)')],
).opts(
    size=1,
    alpha=0.3,
)

points.hist(dimension=['Kd', 'Kdf'])

To visualize the theoretical curve, we can perform a posterior predictive check. We plot the middle 95th percentile of anisotropy curves acquired from theoretical curves arising from our sampled parameters. This requires repeated calls to the anisotropy() function, which in turn calls eqtk.solveNK(). First, we generate smooth sets of initial concentrations to use in generating the plot.

c0_d_smooth = np.zeros((100, 5))
c0_d_smooth[:, 0] = np.linspace(0, 35, 100)
c0_d_smooth[:, 2] = c0_d[0, 2]
x0_d_smooth = c0_d_smooth / solvent_density

c0_c_smooth = np.zeros((100, 5))
c0_c_smooth[:, 0] = c0_c[0, 0]
c0_c_smooth[:, 1] = np.linspace(0, 200, 100)
c0_c_smooth[:, 2] = c0_c[0, 2]
x0_c_smooth = c0_c_smooth / solvent_density

Next, we extract parameter values from the thinned samples.

logKd = df_samples['logKd'].values[::10]
logKdf = df_samples['logKdf'].values[::10]
rf_d = df_samples['rf_d'].values[::10]
rf_c = df_samples['rf_c'].values[::10]
rb_d = df_samples['rb_d'].values[::10]
rb_c = df_samples['rb_c'].values[::10]

And finally, we can compute the posterior predictive distribution of the dissociation and competition anisotropy curves.

dissociation_ppc = np.empty((len(logKd), len(x0_d_smooth)))
competition_ppc = np.empty((len(logKd), len(x0_c_smooth)))
for i in range(len(logKd)):
    dissociation_ppc[i, :] = anisotropy(
        x0_d_smooth, N, logKd[i], logKdf[i], rf_d[i], rb_d[i]
    )
    competition_ppc[i, :] = anisotropy(
        x0_c_smooth, N, logKd[i], logKdf[i], rf_c[i], rb_c[i]
    )

With the curves in hand, we can compute the relevant percentiles and make the plots.

# Compute percentiles from curves
r_d_ppc = np.percentile(dissociation_ppc, [2.5, 50, 97.5], axis=0)
r_c_ppc = np.percentile(competition_ppc, [2.5, 50, 97.5], axis=0)

# Fluorophore dissociation curve
qtile_60_d = hv.Area(
    data=(c0_d_smooth[:, 0], r_d_ppc[2], r_d_ppc[0]),
    kdims="[actin] (µM)",
    vdims=["anisotropy", "dummy"],
    label="fluorophore dissociation",
).opts(line_color=None, fill_color="orange", fill_alpha=0.3)

median_curve_d = hv.Curve(
    data=(c0_d_smooth[:, 0], r_d_ppc[1]), label="fluorophore dissociation"
).opts(color="orange")

diss_ppc_plot = (qtile_60_d * median_curve_d * diss_plot).opts(
    show_legend=False
)

# Competition curve curve
qtile_60_c = hv.Area(
    data=(c0_c_smooth[:, 1], r_c_ppc[2], r_c_ppc[0]),
    kdims="[SD] (µM)",
    vdims=["anisotropy", "dummy"],
    label="competition",
).opts(line_color=None, fill_color="orange", fill_alpha=0.3)

median_curve_c = hv.Curve(
    data=(c0_c_smooth[:, 1], r_c_ppc[1]), label="competition"
).opts(color="orange")

comp_ppc_plot = (qtile_60_c * median_curve_c * comp_plot).opts(
    show_legend=False
)

# Construct layout
(diss_ppc_plot + comp_ppc_plot).opts(shared_axes=False)

Metropolis sampling algorithm

Below is the source code for the JITted Metropolis sampling used in this case study.

import numpy as np
import pandas as pd
import numba


@numba.jit(nopython=True)
def _adjust_sigma(acc_rate, sigma):
    """
    Tune sigma in proposal distribution.

    Parameters
    ----------
    acc_rate : float
        The acceptance rate.
    sigma : ndarray
        Array of standard deviations for Gaussian proposal
        distribution.

    Returns
    -------
    output : ndarray
        Updated `sigma` values.
    """
    if acc_rate < 0.001:
        return sigma * 0.1
    elif acc_rate < 0.05:
        return sigma * 0.5
    elif acc_rate < 0.2:
        return sigma * 0.9
    elif acc_rate > 0.95:
        return sigma * 10.0
    elif acc_rate > 0.75:
        return sigma * 2.0
    elif acc_rate > 0.5:
        return sigma * 1.1
    else:
        return sigma


def mh_sample(logpost, x0, sigma=None, discrete=None, args=(), n_burn=1000,
              n_steps=1000, tune_interval=100, variable_names=None,
              return_acceptance_rate=False):
    """
    Parameters
    ----------
    logpost : function
        The function to compute the log posterior. It has call
        signature `logpost(x, *args)`, and must by compiled with
        numba in nopython mode.
    x0 : ndarray, shape (n_variables,)
        The starting location of a walker in parameter space.
    sigma : ndarray, shape (n_variables,)
        The standard deviations for the proposal distribution.
        If None, takes all values to be one.
    discrete : ndarray of bools, shape (n_variables,)
        discrete[i] is True if variable i is discrete and False
        otherwise. If None (default), all entries are False.
    args : tuple
        Additional arguments passed to `logpost()` function.
    n_burn : int, default 1000
        Number of burn-in steps.
    n_steps : int, default 1000
        Number of steps to take after burn-in.
    tune_interval : int, default 100
        Number of steps to use when determining acceptance
        fraction for tuning.
    variable_names : list, length n_variables
        List of names of variables. If None, then variable names
        are sequential integers.
    return_acceptance_rate : bool, default False
        If True, also return acceptance rate.

    Returns
    -------
    output : DataFrame
        The first `n_variables` columns contain the samples.
        Additionally, column 'lnprob' has the log posterior value
        at each sample.
    """

    @numba.jit(nopython=True)
    def mh_step(x, logpost_current, sigma, discrete, args=()):
        """
        Parameters
        ----------
        x : ndarray, shape (n_variables,)
            The present location of the walker in parameter space.
        logpost_current : float
            The current value of the log posterior.
        sigma : ndarray, shape (n_variables, )
            The standard deviations for the proposal distribution.
        discrete : ndarray of bools, shape (n_variables,)
            discrete[i] is True if variable i is discrete and False
            otherwise. If None (default), all entries are False.
        args : tuple
            Additional arguments passed to `logpost()` function.

        Returns
        -------
        output : ndarray, shape (n_variables,)
            The position of the walker after the Metropolis-Hastings
            step. If no step is taken, returns the inputted `x`.
        """
        # Draw the next step
        x_next = np.empty(len(x), dtype=x.dtype)
        for i, x_val in enumerate(x):
            if discrete[i]:
                x_next[i] = np.round(np.random.normal(x[i], sigma[i]))
            else:
                x_next[i] = np.random.normal(x[i], sigma[i])

        # Compute log posterior
        logpost_new = logpost(x_next, *args)

        # Compute the log Metropolis ratio
        log_r = logpost_new - logpost_current

        # Accept or reject step
        if log_r >= 0 or np.random.random() < np.exp(log_r):
            return x_next, logpost_new, True
        else:
            return x, logpost_current, False

        # Never tune if asked not to
        if tune_interval <= 0:
            tune_interval = n_burn + n_steps
        elif tune_interval < 50:
            raise RuntimeError('Tune interval should be at least 50.')


    @numba.jit(nopython=True)
    def _metropolis_sample(x0, sigma, discrete, n_burn, n_steps,
                           tune_interval=100, args=()):
        """
        Numba'd sampler.
        """

        # Initialize
        x = np.copy(x0)
        n_accept = 0
        n_accept_total = 0
        n = 0
        n_tune_steps = 0
        n_continuous_steps = 0
        lnprob = np.empty(n_steps)
        logpost_current = logpost(x, *args)

        # Burn in
        while n < n_burn:
            while n_tune_steps < tune_interval and n < n_burn:
                x, logpost_current, accept = mh_step(
                            x, logpost_current, sigma, discrete, args)
                n += 1
                n_tune_steps += 1
                n_accept += accept
            sigma = _adjust_sigma(n_accept/tune_interval, sigma)
            n_accept = 0
            n_tune_steps = 0

        # Samples
        x_samples = np.empty((n_steps, len(x)))
        n = 0
        while n < n_steps:
            while n_tune_steps < tune_interval and n < n_steps:
                x, logpost_current, accept = mh_step(
                            x, logpost_current, sigma, discrete, args)
                x_samples[n,:] = x
                lnprob[n] = logpost_current
                n_tune_steps += 1
                n_accept += accept
                n_accept_total += accept
                n += 1
            sigma = _adjust_sigma(n_accept/tune_interval, sigma)
            n_accept = 0
            n_tune_steps = 0

        return x_samples, lnprob, n_accept_total/n_steps

    # Use default indices for variable names
    if variable_names is None:
        variable_names = list(range(len(x0)))

    # Set sigma to be default if not provided
    if sigma is None:
        sigma = np.ones(len(x0))

    # Default to continuous variables
    if discrete is None:
        discrete = np.array([False]*len(x0))
    else:
        discrete = np.array(discrete)

    # Grab the samples in a NumPy array
    x_samples, lnprob, acc_rate = _metropolis_sample(
            x0, sigma, discrete, n_burn, n_steps,
            tune_interval=tune_interval, args=args)

    df = pd.DataFrame(columns=variable_names, data=x_samples)
    df['lp__'] = lnprob

    if return_acceptance_rate:
        return df, acc_rate
    return df