fit-binding-models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "70836125-de12-410d-8b6b-2d883b74596b",
   "metadata": {},
   "source": [
    "## Fitting models to data\n",
    "\n",
    "In biochemical studies, we often want to infer quantiative information about reactions from experiments. This often requires fitting a quantiative model to data. In this notebook, we are going to explore how this process works, how to assess the results, and what can go wrong. \n",
    "\n",
    "### The reaction\n",
    "\n",
    "Consider a binding reaction between a macromolecule $M$ and a small molecule $X$:\n",
    "\n",
    "$$MX \\rightleftarrows M + X$$\n",
    "\n",
    "We are often interested in the dissociation constant ($K_{D}$). This is defined as:\n",
    "\n",
    "$$K_{D} \\equiv \\frac{[M][X]}{[MX]}.$$\n",
    "\n",
    "To estimate $K_{D}$, one often measures some experimental observable as a function of the concentration of $[X]$ added to solution. In a perfect experiment, one would directly measure the fraction of macromolecule bound to $X$:\n",
    "\n",
    "$$\\Theta = \\frac{[MX]}{[M]+[MX]}.$$\n",
    "\n",
    "### Developing a mathematical model\n",
    "\n",
    "A single-site binding reaction can be described with the following simple mathematical model:\n",
    "\n",
    "$$\\Theta \\sim \\frac{[X]/K_{D}}{1 + [X]/K_{D}}.$$\n",
    "\n",
    "One can rarely measure $\\Theta$ directly. Instead, we often have some observable that changes when $X$ binds to $M$. Consider an extreme case where a protein is unfolded when unbound, then folds when bound. We could use circular dichroism to measure the formation of secondary structure as a function of $X$, then relate this to $\\Theta$. The function would look something like the following. Our observed signal at some concentration of $X$ ($S_{obs}$) arises from a mixture of the signals from $M$ ($S_{M}$) and $MX$ ($S_{MX}$). \n",
    "\n",
    "$$S_{obs} \\sim S_{M}(1-\\Theta) + S_{MX}\\Theta$$\n",
    "\n",
    "Think about this for a moment. If the protein is all bound, $\\Theta$ will be 1.0 and $S_{obs}$ will equal $S_{MX}$. If none of the protein is bound, $\\Theta$ will be 0.0 and $S_{obs}$ will equal $S_{M}$. At intermediate values of $\\Theta$, we will have proportional mixtures of $S_{M}$ and $S_{MX}$. \n",
    "\n",
    "We can rearrange this to:\n",
    "\n",
    "$$S_{obs} \\sim S_{M} + \\Theta (S_{MX} - S_{M})$$\n",
    "\n",
    "We can substitute $\\Theta$ from above, giving our final model:\n",
    "\n",
    "$$S_{obs} \\sim f(K_{D},S_{M},S_{MX},[X]) = S_{M} + \\frac{[X]/K_{D}}{1 + [X]/K_{D}} (S_{MX} - S_{M})$$\n",
    "\n",
    "This function $f$ has three fit parameters ($K_{D}$, $S_{M}$ and $S_{MX}$) and one control variable ($[X]$). The parameter we care about is $K_{D}$; the other two parameters let us map binding to our signal $S_{obs}$. We can experimentally manipulate the control variable $[X]$ and measure $S_{obs}$, thus allowing us to learn about $K_{D}$. \n",
    "\n",
    "### Estimating the parameters of the model\n",
    "\n",
    "We now need to find the parameters of this model that reproduce our observed data. To do so, we use [nonlinear regression](https://en.wikipedia.org/wiki/Nonlinear_regression). This is a glorified form of guess-and-check, where we try to identify values for the parameters that minimize the difference between our model and observed data. In most nonlinear regression, we minimize the *least squares difference* between our calculated and observed data. If we have measured $S_{obs}$ at $N$ different concentrations of $[X]$, we would write this as:\n",
    "\n",
    "$$D = \\sum_{i=1}^{i \\le N} \\Big ( f(K_{D},S_{M},S_{MX},[X]_{i}) - S_{obs,[X]_{i}}) \\Big )^{2}.$$\n",
    "\n",
    "When we minimize $D$ with respect to $K_{D}$, $S_{M}$, and $S_{MX}$, we obtain the *maximum likelihood estimate* of our parameter of interest, $K_{D}$. \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "960064ad-d25c-4efc-9c02-5c8cb42da0b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@title Initialize environment\n",
    "\n",
    "#@markdown We recommend setting up a working directory on your google drive. This is a \n",
    "#@markdown convenient way to pass files in and out of this analysis. It will \n",
    "#@markdown also allow you to save your work. If you put `biophysics` into the form\n",
    "#@markdown field below, the analyis will save all of its calculations in the \n",
    "#@markdown `biophysics` directory in MyDrive (i.e. the top directory at\n",
    "#@markdown https://drive.google.com). This script will create the directory if \n",
    "#@markdown it does not already exist. If the directory already exists, any files\n",
    "#@markdown that are already in that directory will be available for the analysis. \n",
    "#@markdown You could, for example, put a file called `data.csv` in `biophysics` and then\n",
    "#@markdown access it as \"data.csv\" in all cells below.\n",
    "#@markdown <br/>\n",
    "#@markdown Note: Google may prompt you for permission to access the drive. \n",
    "#@markdown To work in a temporary colab environment, leave this blank. Your results\n",
    "#@markdown will disappear when you close the directory. \n",
    "\n",
    "try:\n",
    "    import google.colab\n",
    "    RUNNING_IN_COLAB = True\n",
    "except ImportError:\n",
    "    RUNNING_IN_COLAB = False\n",
    "except Exception as e: \n",
    "    err = \"Could not figure out if runnning in a colab notebook\\n\"\n",
    "    raise Exception(err) from e\n",
    "\n",
    "# ------------------------------------------------------------------------------\n",
    "# Imports\n",
    "\n",
    "if RUNNING_IN_COLAB:\n",
    "    %pip install -q ipywidgets\n",
    "\n",
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy\n",
    "from scipy import optimize\n",
    "\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "\n",
    "import os\n",
    "import colorsys\n",
    "\n",
    "# ------------------------------------------------------------------------------\n",
    "# Environment\n",
    "\n",
    "if RUNNING_IN_COLAB:\n",
    "    \n",
    "    working_dir = \"/content/\"\n",
    "\n",
    "    # Select a working directory on google drive\n",
    "    google_drive_directory = \"\" #@param {type:\"string\"}\n",
    "    google_drive_directory = google_drive_directory.strip()\n",
    "\n",
    "    # Set up google drive\n",
    "    if google_drive_directory != \"\":\n",
    "\n",
    "        from google.colab import drive\n",
    "        drive.mount('/content/gdrive/')\n",
    "\n",
    "        working_dir = f\"/content/gdrive/MyDrive/{google_drive_directory}\"\n",
    "        os.system(f\"mkdir -p {working_dir}\")\n",
    "\n",
    "    os.chdir(working_dir)\n",
    "    print(f\"Working directory: {os.getcwd()}/\")\n",
    "\n",
    "    print(\"\\nCurrent directory contents:\")\n",
    "    print(os.getcwd())\n",
    "    for f in os.listdir(\".\"):\n",
    "        print(f\"    {f}\")\n",
    "    print()\n",
    "    \n",
    "\n",
    "# ------------------------------------------------------------------------------\n",
    "# Default graph label sizing\n",
    "\n",
    "SMALL_SIZE = 14\n",
    "MEDIUM_SIZE = 16\n",
    "BIGGER_SIZE = 18\n",
    "\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "\n",
    "# ------------------------------------------------------------------------------\n",
    "# Fitting functions\n",
    "\n",
    "def binding_model_full(Kd,Sm,Smx,Xtot,Mtot):\n",
    "    \"\"\"\n",
    "    Single site binding model (SM <--> S + M) with no assumption that Mtot is\n",
    "    less than Kd. \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    Kd : float\n",
    "        dissociation constant\n",
    "    Sm : float\n",
    "        signal of the M form \n",
    "    Smx : float\n",
    "        signal of the MX form\n",
    "    Xtot : float\n",
    "        total X concentration\n",
    "    Mtot : float\n",
    "        total M concentration\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    signal : float\n",
    "        spectroscopic signal given the current [MX]/([M] + [MX])\n",
    "    \"\"\"\n",
    "    \n",
    "    a = 1\n",
    "    b = -(Xtot + Mtot + Kd)\n",
    "    c = Mtot*Xtot\n",
    "    \n",
    "    MX = (-b - np.sqrt(b**2 - 4*a*c))/(2*a)\n",
    "    \n",
    "    theta = MX/Mtot\n",
    "    \n",
    "    signal = Sm + theta*(Smx - Sm)\n",
    "    \n",
    "    return signal\n",
    "\n",
    "def binding_model(Kd,Sm,Smx,Xtot):\n",
    "    \"\"\"\n",
    "    Single site binding model (SM <--> S + M)\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    Kd : float\n",
    "        dissociation constant\n",
    "    Sm : float\n",
    "        signal of the M form \n",
    "    Smx : float\n",
    "        signal of the MX form\n",
    "    X : float\n",
    "        total X concentration\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    signal : float\n",
    "        spectroscopic signal given the current [MX]/([M] + [MX])\n",
    "    \"\"\"\n",
    "    \n",
    "    alpha = Xtot/Kd\n",
    "    theta = alpha/(1 + alpha)\n",
    "    \n",
    "    obs = Sm + theta*(Smx - Sm)\n",
    "    \n",
    "    return obs\n",
    "\n",
    "def binding_model_residuals(param,x_expt,y_expt,Mtot=None):\n",
    "    \"\"\"\n",
    "    Calculate the difference between the signal calculated from the model and\n",
    "    the observed data.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    param : list-like of floats\n",
    "        model parameters (Kd, Sm, Smx)\n",
    "    x_expt : list-like of floats\n",
    "        Xtot values at which data were collected\n",
    "    y_expt : list-like of floats\n",
    "        observed signal at the Xtot values given in x_expt\n",
    "    Mtot : float, optional\n",
    "        total M concentration. (If specified, use full binding model; if not,\n",
    "        used reduced model that assumes [M]tot << Kd\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    residual : np.ndarray\n",
    "        difference between calculated signal and y_expt for each value of \n",
    "        x_expt\n",
    "    \"\"\"\n",
    "    \n",
    "    x_expt = np.array(x_expt)\n",
    "    y_expt = np.array(y_expt)\n",
    "    \n",
    "    if Mtot is None:\n",
    "        calc = binding_model(param[0],param[1],param[2],x_expt)\n",
    "    else:\n",
    "        calc = binding_model(param[0],param[1],param[2],x_expt,Mtot)\n",
    "    \n",
    "    residual = calc - y_expt\n",
    "    \n",
    "    return residual\n",
    "\n",
    "\n",
    "def _get_R2(y_calc,y_expt):\n",
    "    \n",
    "    RSS = np.sum((y_expt - y_calc)**2)\n",
    "    TSS = np.sum((y_expt - np.mean(y_expt))**2)\n",
    "    \n",
    "    return 1 - RSS/TSS\n",
    "\n",
    "def _run_regresssion(fcn,\n",
    "                     residual_fcn,\n",
    "                     param_guesses,\n",
    "                     x_expt,\n",
    "                     y_expt,\n",
    "                     Mtot=None,\n",
    "                     ax=None,\n",
    "                     residual_ax=None,\n",
    "                     color=\"black\",\n",
    "                     alpha=1.0,\n",
    "                     max_nfev=None,\n",
    "                     label=None):\n",
    "    \"\"\"\n",
    "    Use nonlinear regression to fit a model to data.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    fcn : function\n",
    "        function to analyze\n",
    "    residual_fcn : function\n",
    "        residual function\n",
    "    param_guesses : list-like of floats\n",
    "        model parameters\n",
    "    x_expt : list-like of floats\n",
    "        Xtot values at which data were collected\n",
    "    y_expt : list-like of floats\n",
    "        observed signal at the Xtot values given in x_expt\n",
    "    Mtot : float, optional\n",
    "        total M concentration. (If specified, use full binding model; if not,\n",
    "        used reduced model that assumes [M]tot << Kd\n",
    "    ax : matplotlib.Axis, optional\n",
    "        if specified, draw the points and the fit on this ax\n",
    "    residual_ax : matplotlib.Axis, optional\n",
    "        if specified, draw the residual points on this axis\n",
    "    color : str, default=\"black\"\n",
    "        color to use to draw series\n",
    "    alpha : float, default=1.0\n",
    "        opacity to give to series\n",
    "    max_nfev : int, optional\n",
    "        maximum number of iterations to run the regression\n",
    "    label : str, optional\n",
    "        label to assign to the newly drawn series\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    fit.x : np.ndarray\n",
    "        maximum likelihood fit parameter estimates\n",
    "    x_calc : np.ndarray\n",
    "        values of x at which calculation is done to get smooth line\n",
    "    y_calc : np.ndarray\n",
    "        function values calculated at x_calc\n",
    "    y_calc_expt : np.ndarray\n",
    "        function values calcualated at x_expt\n",
    "    \"\"\"\n",
    "    \n",
    "    x_min = np.min(x_expt)\n",
    "    x_max = np.max(x_expt)\n",
    "    y_min = np.min(y_expt)\n",
    "    y_max = np.max(y_expt)\n",
    "    del_y = y_max - y_min\n",
    "    x_calc = np.linspace(x_min*.9,x_max*1.1,100)\n",
    "    \n",
    "    fit = optimize.least_squares(residual_fcn,\n",
    "                                 x0=param_guesses,\n",
    "                                 kwargs={\"x_expt\":x_expt,\n",
    "                                         \"y_expt\":y_expt,\n",
    "                                         \"Mtot\":Mtot},\n",
    "                                 max_nfev=max_nfev)\n",
    "    \n",
    "    y_calc = fcn(*fit.x,x_calc)\n",
    "    y_calc_expt = fcn(*fit.x,x_expt)\n",
    "            \n",
    "    if ax is not None:\n",
    "        ax.plot(x_calc,y_calc,'-',lw=2,color=color,alpha=alpha,label=label)\n",
    "\n",
    "    if residual_ax is not None:\n",
    "\n",
    "        residual_ax.plot(x_expt,y_calc_expt - y_expt,'o',color=color,alpha=alpha)\n",
    "            \n",
    "    return fit.x, x_calc, y_calc, y_calc_expt\n",
    "    \n",
    "def fit_data(x_expt,\n",
    "             y_expt,\n",
    "             y_err,\n",
    "             Kd_guess,\n",
    "             Sm_guess,\n",
    "             Smx_guess,\n",
    "             Mtot=None,\n",
    "             num_bootstraps=0):\n",
    "    \"\"\"\n",
    "    Fit a binding model to the data given in x_expt, y_expt, and y_err. Generate\n",
    "    plots and potentially bootstrap to assess parameter uncertainty. Model has \n",
    "    the form: MX <--> M + X. We assume there is a signal Sm proportional to [M]\n",
    "    and another signal Smx proportional to [MX]. \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    x_expt : list-like of floats\n",
    "        Xtot values at which data were collected\n",
    "    y_expt : list-like of floats\n",
    "        observed signal at the Xtot values given in x_expt\n",
    "    y_err : list-like of floats\n",
    "        standard deviation of measurements for each y_expt\n",
    "    Kd_guess : float\n",
    "        guess for the Kd\n",
    "    Sm_guess : float\n",
    "        guess for the signal of M\n",
    "    Smx_guess : float\n",
    "        guess for the signal of MX\n",
    "    Mtot : float, optional\n",
    "        total M concentration. (If specified, use full binding model; if not,\n",
    "        used reduced model that assumes [M]tot << Kd \n",
    "    num_bootstraps : int, default=0\n",
    "        do this many bootstrap pseudreplicate fits\n",
    "    \"\"\"\n",
    "            \n",
    "    # Get data limits\n",
    "    x_min = np.min(x_expt)\n",
    "    x_max = np.max(x_expt)\n",
    "    y_min = np.min(y_expt)\n",
    "    y_max = np.max(y_expt)\n",
    "    del_y = y_max - y_min\n",
    "        \n",
    "    # Create plot\n",
    "    fig, ax = plt.subplots(2,1,figsize=(6,12))\n",
    "    \n",
    "    # Create data plot\n",
    "    ax[0].plot(x_expt,y_expt,'o',color=\"black\")\n",
    "    ax[0].errorbar(x_expt,y_expt,y_err,fmt='o',capsize=5,ms=0,lw=1,color=\"black\")\n",
    "    ax[0].set_xlabel(\"$[X]_{tot}$ ($\\mu M$)\")\n",
    "    ax[0].set_ylabel(\"signal\")\n",
    "    ax[0].spines['top'].set_visible(False)\n",
    "    ax[0].spines['right'].set_visible(False)\n",
    "    \n",
    "    # Create residuals plot\n",
    "    ax[1].plot((x_min,x_max),(0,0),'--',lw=2,color=\"gray\")\n",
    "    ax[1].set_xlabel(\"$[X]_{tot}$ ($\\mu M$)\")\n",
    "    ax[1].set_ylabel(\"calc - obs\")\n",
    "    ax[1].spines['top'].set_visible(False)\n",
    "    ax[1].spines['right'].set_visible(False)\n",
    "\n",
    "    # Figure out which model to use\n",
    "    if Mtot is None:\n",
    "        fcn = binding_model\n",
    "    else:\n",
    "        fcn = binding_model_full\n",
    "    \n",
    "    # Do maximum likelhood fit\n",
    "    fit_param, x_calc, y_calc, y_calc_expt = _run_regresssion(fcn,\n",
    "                                                              binding_model_residuals,\n",
    "                                                              [Kd_guess,Sm_guess,Smx_guess],\n",
    "                                                              x_expt,\n",
    "                                                              y_expt,\n",
    "                                                              Mtot,\n",
    "                                                              ax=ax[0],\n",
    "                                                              residual_ax=ax[1],\n",
    "                                                              color=\"black\")\n",
    "    \n",
    "    # Get fit parameters\n",
    "    Kd = fit_param[0]\n",
    "    Sm = fit_param[1]\n",
    "    Smx = fit_param[2]\n",
    "    \n",
    "    R2 = _get_R2(y_calc_expt,y_expt)\n",
    "    \n",
    "    # Figure out y limits on residuals plot\n",
    "    biggest_diff = 1.1*np.max(np.abs(y_calc_expt - y_expt))\n",
    "    \n",
    "    if biggest_diff < del_y/4:\n",
    "        span = del_y/4\n",
    "    else:\n",
    "        span = biggest_diff\n",
    "    \n",
    "    ax[1].set_ylim(-span,span)\n",
    "\n",
    "    \n",
    "    Kd_err = None\n",
    "    if num_bootstraps > 0:\n",
    "    \n",
    "        Kd_err = []\n",
    "        for i in range(num_bootstraps):\n",
    "            \n",
    "            # Generate a pseudoreplicate\n",
    "            this_y_expt = y_expt + np.random.normal(0,y_err,len(y_expt))\n",
    "\n",
    "            # Generate color\n",
    "            s = i/(num_bootstraps)*.9 + 0.1\n",
    "            color = [s,s,1]\n",
    "            \n",
    "            # Plot pseudoreplicate\n",
    "            ax[0].plot(x_expt,this_y_expt,'o',color=color,alpha=0.5)\n",
    "            \n",
    "            # Do pseudoreplicate fit\n",
    "            fit_param, x_calc, y_calc, y_calc_expt = _run_regresssion(binding_model,\n",
    "                                                                      binding_model_residuals,\n",
    "                                                                      [Kd_guess,Sm_guess,Smx_guess],\n",
    "                                                                      x_expt,\n",
    "                                                                      this_y_expt,\n",
    "                                                                      ax=ax[0],\n",
    "                                                                      residual_ax=ax[1],\n",
    "                                                                      color=color,\n",
    "                                                                      alpha=0.5)\n",
    "            \n",
    "            # Record pseudoreplicate Kd\n",
    "            Kd_err.append(fit_param[0])\n",
    "            \n",
    "        # Record error bit\n",
    "        Kd_err = np.std(Kd_err)\n",
    "        Kd_err = f\" $\\pm$ {Kd_err:.1f}\"\n",
    "        \n",
    "    if Kd_err is None:\n",
    "        Kd_err = \"\"\n",
    "        \n",
    "    ax[0].text(x_max*0.85,y_min + del_y*0.15,f\"$K_D$: {Kd:.2f}{Kd_err} $\\mu M$\")\n",
    "    ax[0].text(x_max*0.85,y_min + del_y*0.05,f\"$R^2$: {R2:.4f}\")\n",
    "    \n",
    "# ------------------------------------------------------------------------------\n",
    "# Datasets\n",
    "\n",
    "def read_dataset(origin):\n",
    "    df = pd.read_csv(origin)\n",
    "    try:\n",
    "        df = df.drop(\"Unnamed: 0\",axis=1)\n",
    "    except KeyError:\n",
    "        pass\n",
    "\n",
    "    return df\n",
    "\n",
    "df1 = read_dataset(\"https://raw.githubusercontent.com/harmsm/biochem-jupyter-notebooks/master/binding/raw-binding-data.csv\")\n",
    "df2 = read_dataset(\"https://raw.githubusercontent.com/harmsm/biochem-jupyter-notebooks/master/binding/raw-binding-data_1.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c464cd4f-6759-4d61-a3e6-befdd07a68b7",
   "metadata": {},
   "source": [
    "## What is fitting?\n",
    "\n",
    "The following cell lets you fit the model above to some toy data, with increasing numbers of steps in the regression. (Think, an increasing number of guess-and-check steps). What happens to the parameters as more steps occur? When do you think we should stop doing regression steps? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3b0bc7b-43da-4d94-8709-c1ab11f2127c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@title Press the \"Play\" button on the left to run.\n",
    "\n",
    "def plot_fitting_vs_steps(max_steps):\n",
    "    \n",
    "    Kd_guess = 5\n",
    "    Sm_guess = 0\n",
    "    Smx_guess = 2\n",
    "    x_expt = np.arange(5)\n",
    "    y_expt = np.array([0,1,1.5,1.6,1.7])\n",
    "    y_err = 0.1*np.ones(5)\n",
    "\n",
    "    fig, ax = plt.subplots(1,figsize=(6,6))\n",
    "\n",
    "    # Create data plot\n",
    "    ax.plot(x_expt,y_expt,'o',color=\"black\",zorder=100)\n",
    "    ax.errorbar(x_expt,y_expt,y_err,fmt='o',capsize=5,ms=0,lw=1,color=\"black\",zorder=100)\n",
    "    ax.set_xlabel(\"$[X]_{tot}$ ($\\mu M$)\")\n",
    "    ax.set_ylabel(\"signal\")\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.set_ylim([0,2])\n",
    "    \n",
    "    # Do maximum likelhood fit\n",
    "    for i in range(1,max_steps+1):\n",
    "\n",
    "        color = colorsys.hsv_to_rgb(i/10,1,0.8)\n",
    "\n",
    "        fit_param, x_calc, y_calc, y_calc_expt = _run_regresssion(binding_model,\n",
    "                                                                  binding_model_residuals,\n",
    "                                                                  [Kd_guess,Sm_guess,Smx_guess],\n",
    "                                                                  x_expt,\n",
    "                                                                  y_expt,\n",
    "                                                                  ax=ax,\n",
    "                                                                  color=color,\n",
    "                                                                  max_nfev=i+1,\n",
    "                                                                  label=f\"step {i}\")\n",
    "        print(f\"Step: {i}, Kd: {fit_param[0]:.2f}, Sm: {fit_param[1]:.2f}, Smx: {fit_param[2]:.2f}\")\n",
    "\n",
    "\n",
    "    ax.legend(loc=\"lower right\")\n",
    "    \n",
    "#Kd_guess_slider = widgets.IntSlider(min=1,max=10,by=1,value=5,description=\"guess Kd (micromolar)\")\n",
    "max_steps_slider = widgets.IntSlider(min=1,max=10,by=1,value=1,description=\"max_steps\")\n",
    "\n",
    "w = widgets.interactive(plot_fitting_vs_steps,\n",
    "                        #Kd_guess=Kd_guess_slider,\n",
    "                        max_steps=max_steps_slider)\n",
    "                        \n",
    "display(w)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb0ebcdd-2262-42ab-b0ef-109921877f04",
   "metadata": {},
   "source": [
    "## Initial parameter guesses\n",
    "\n",
    "The result of a fit may depend on our initial guesses for our fit parameter values. A good rule of thumb is to choose plausible first guess, as it makes it more likely you find the relevant solution. Can you make the results change by changing your parameter guesses? \n",
    "\n",
    "(Does this rule of thumb worry you? It should... Why?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb0d19a1-8981-4f9c-aa3b-5733f7247115",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@title Press the \"Play\" button on the left to run.\n",
    "\n",
    "Kd_guess_slider = widgets.FloatSlider(min=-100,max=100,by=5,value=10,description=\"guess Kd (micromolar)\")\n",
    "Sm_guess_slider = widgets.FloatSlider(min=-100,max=100,by=5,value=0,description=\"guess Sm\")\n",
    "Smx_guess_slider = widgets.FloatSlider(min=-100,max=100,by=5,value=2,description=\"guess Smx (micromolar)\")\n",
    "\n",
    "def _fit_data_wrapper(Kd_guess,Sm_guess,Smx_guess):\n",
    "    \n",
    "    return fit_data(x_expt=np.arange(5),\n",
    "                    y_expt=np.array([0,1,1.5,1.6,1.7]),\n",
    "                    y_err=0.1*np.ones(5),\n",
    "                    Kd_guess=Kd_guess,\n",
    "                    Sm_guess=Sm_guess,\n",
    "                    Smx_guess=Smx_guess,\n",
    "                    num_bootstraps=0)\n",
    "    \n",
    "\n",
    "w = widgets.interactive(_fit_data_wrapper,\n",
    "                        Kd_guess=Kd_guess_slider,\n",
    "                        Sm_guess=Sm_guess_slider,\n",
    "                        Smx_guess=Smx_guess_slider)\n",
    "                        \n",
    "display(w)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ef12b4a-da32-40a3-a88a-dff896be6cc4",
   "metadata": {},
   "source": [
    "## How do we assess fit quality?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4aa2764-2029-4c21-8474-5744cbce2a5d",
   "metadata": {},
   "source": [
    "How do we decide if our fit is good? There are two major metrics: $R^{2}$ and the *fit residuals*. \n",
    "\n",
    "+ $R^{2}$ measures how much scatter there is off the line. If the line drawn by our model goes through every point exactly, $R^{2}$ would be 1: the fit is perfect. If the line is far away from the points, $R^{2}$ will be close to zero. This value is shown on the top plot. \n",
    "+ The *fit residuals* determine if the model goes through the data with the same level of \"goodness\" as a function of our control variable ($[X]$). This is the middle plot. It is calculated as $f([X]_{i}) - S_{obs,i}$. A good fit will have *random* residuals, with a similar scatter across all values of $[X]$. If you make a histogram of these values (bottom plot), it should be normally distributed and centered at zero. \n",
    "\n",
    "Compare the values of $R^{2}$ and the residuals using either the \"Good data\" or not (click toggle). Why is the fit to the good data more valid than the fit to the not-so-good data?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "383c425f-9dd3-4192-b222-1c605a97e5dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@title Press the \"Play\" button on the left to run.\n",
    "\n",
    "\n",
    "def fit_quality(good_model=True):\n",
    "    \n",
    "    Kd_guess = 5\n",
    "    Sm_guess = 0\n",
    "    Smx_guess = 2\n",
    "    \n",
    "    x_expt = np.linspace(0,30,61)    \n",
    "    if good_model:\n",
    "        y_expt = binding_model(5,0,2,x_expt) + np.random.normal(0,0.1,len(x_expt))\n",
    "        y_err = 0.1*np.ones(len(x_expt))\n",
    "    else:\n",
    "        y_expt = x_expt**1.3 + np.random.normal(0,0.1,len(x_expt))\n",
    "        y_err = 1*np.ones(len(x_expt))\n",
    "    \n",
    "\n",
    "    # Create plot\n",
    "    fig, ax = plt.subplots(3,1,figsize=(6,18))\n",
    "    \n",
    "    # Create data plot\n",
    "    ax[0].plot(x_expt,y_expt,'o',color=\"black\")\n",
    "    ax[0].errorbar(x_expt,y_expt,y_err,fmt='o',capsize=5,ms=0,lw=1,color=\"black\")\n",
    "    ax[0].set_xlabel(\"$[X]_{tot}$ ($\\mu M$)\")\n",
    "    ax[0].set_ylabel(\"signal\")\n",
    "    ax[0].spines['top'].set_visible(False)\n",
    "    ax[0].spines['right'].set_visible(False)\n",
    "        \n",
    "        \n",
    "    # Create residuals plot\n",
    "    ax[1].plot((0,30*1.1),(0,0),'--',lw=2,color=\"gray\")\n",
    "    ax[1].set_xlabel(\"$[X]_{tot}$ ($\\mu M$)\")\n",
    "    ax[1].set_ylabel(\"calc - obs\")\n",
    "    ax[1].spines['top'].set_visible(False)\n",
    "    ax[1].spines['right'].set_visible(False)\n",
    "    \n",
    "    fit_param, x_calc, y_calc, y_calc_expt = _run_regresssion(binding_model,\n",
    "                                                              binding_model_residuals,\n",
    "                                                              [Kd_guess,Sm_guess,Smx_guess],\n",
    "                                                              x_expt,\n",
    "                                                              y_expt,\n",
    "                                                              ax=ax[0],\n",
    "                                                              residual_ax=ax[1])\n",
    "    \n",
    "    if good_model:\n",
    "        ax[2].hist(y_calc_expt-y_expt,bins=np.arange(-0.3,0.38,0.08))\n",
    "    else:\n",
    "        ax[2].hist(y_calc_expt-y_expt)\n",
    "\n",
    "    ax[0].text(25,1,f\"R2: {_get_R2(y_calc_expt,y_expt):.4f}\")\n",
    "\n",
    "\n",
    "good_model_toggle = widgets.ToggleButton(\n",
    "    value=True,\n",
    "    description='Use good data',\n",
    "    disabled=False,\n",
    "    button_style='', # 'success', 'info', 'warning', 'danger' or ''\n",
    "    tooltip='If set to True, use data matching model; otherwise, use non-matching data',\n",
    "    icon='check'\n",
    ")\n",
    "    \n",
    "w = widgets.interactive(fit_quality,\n",
    "                        good_model=good_model_toggle)\n",
    "\n",
    "display(w)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c0abffe-7ba7-4ae6-ab2b-6945a4e9f21b",
   "metadata": {},
   "source": [
    "## How to we assess uncertainty?\n",
    "\n",
    "How do we get error bars on our fit results? One way to do so is using pseudoreplicate bootstrapping. Can you figure out how bootstrapping works playing with the slider below? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57391c1b-4975-4c32-b956-6f39de5c4605",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@title Press the \"Play\" button on the left to run.      \n",
    "\n",
    "num_bootstraps_slider = widgets.IntSlider(min=0,max=10,by=1,value=0,description=\"num bootstraps\")\n",
    "\n",
    "def _fit_data_wrapper(num_bootstraps):\n",
    "    \n",
    "    return fit_data(x_expt=np.arange(5),\n",
    "                    y_expt=np.array([0,1,1.5,1.6,1.7]),\n",
    "                    y_err=0.1*np.ones(5),\n",
    "                    Kd_guess=5,\n",
    "                    Sm_guess=0,\n",
    "                    Smx_guess=2,\n",
    "                    num_bootstraps=num_bootstraps)\n",
    "    \n",
    "\n",
    "w = widgets.interactive(_fit_data_wrapper,\n",
    "                        num_bootstraps=num_bootstraps_slider)\n",
    "                        \n",
    "display(w)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77b71c64-c7da-40bf-9d38-a56aee61812a",
   "metadata": {},
   "source": [
    "## A real example\n",
    "\n",
    "The following two cells show fits to experimental data collected for the binding of two different proteins to a small molecule. The researchers who published these data argue that the protein in the lower panel has a lower affinity for $X$ than the protein in the top panel. Do you believe them? Please justify your answer. What might you do next to further (or correct) their analysis? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9e8b9fe-4a3b-4d68-802c-b1fd3d848890",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@title Press the \"Play\" button on the left to run.   \n",
    "fit_data(df1.Xtot,df1.Sobs,df1.Sobs_err,100,0,1,num_bootstraps=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42ccb75b-6b40-409b-8491-fb24402399fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "#@title Press the \"Play\" button on the left to run.   \n",
    "fit_data(df2.Xtot,df2.Sobs,df2.Sobs_err,100,0,1,num_bootstraps=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "414560bf-882e-4228-8265-6504c30deeae",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}