parameter inference for kinetic model using HMC

induced_fit_parameter_inference.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import theano\n",
    "import theano.tensor as tt\n",
    "import pymc3 as pm\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward_model_numpy(rates, amplitudes, concentrations, time):\n",
    "    L0, P0 = concentrations\n",
    "    k_on, k_off, k_conf_plus, k_conf_minus = rates\n",
    "    c0 = [P0, 0.0, 0.0]\n",
    "    M = np.array([\n",
    "        [-k_on*L0,  k_off, 0.0],\n",
    "        [ k_on*L0, -k_off - k_conf_plus, k_conf_minus],\n",
    "        [ 0.0,      k_conf_plus,        -k_conf_minus]\n",
    "    ])\n",
    "    eigenvalues, V = np.linalg.eig(M)\n",
    "    z = np.linalg.inv(V).dot(c0)\n",
    "    populations = np.einsum('ij,tj,j->ti', V, np.exp(time[:, np.newaxis]*eigenvalues), z)\n",
    "    return populations.dot(amplitudes)\n",
    "\n",
    "def forward_model_theano(rates_, amplitudes_, concentrations_, time):\n",
    "    theano.config.compute_test_value = 'ignore'\n",
    "    concentrations = theano.tensor.dvector('concentrations')\n",
    "    rates = theano.tensor.dvector('rates')\n",
    "    amplitudes = theano.tensor.dvector('amplitudes')\n",
    "    c0 = tt.stacklists([concentrations[1], 0.0, 0.0])\n",
    "    M = tt.stacklists([\n",
    "        [-rates[0]*concentrations[0],  rates[1], 0.0],\n",
    "        [ rates[0]*concentrations[0], -rates[1] - rates[2], rates[3]],\n",
    "        [ 0.0,                         rates[2],           -rates[3]]\n",
    "    ])\n",
    "    eigenvalues, V = tt.nlinalg.eig(M)\n",
    "    V_inv = tt.nlinalg.matrix_inverse(V)\n",
    "    z = tt.dot(V_inv, c0)\n",
    "    a = tt.dot(V.T, amplitudes)\n",
    "    ifmodel = tt.sum(tt.exp(eigenvalues*time[:, np.newaxis])*a*z, axis=1)\n",
    "    ifmodel_function = theano.function([rates, amplitudes, concentrations], ifmodel)\n",
    "    return ifmodel_function(rates_, amplitudes_, concentrations_)\n",
    "\n",
    "k_on = 1.3\n",
    "k_off = 23.0\n",
    "k_conf_plus = 1.3\n",
    "k_conf_minus = 6E-4\n",
    "rates_ = [k_on, k_off, k_conf_plus, k_conf_minus]\n",
    "P0 = 0.1\n",
    "L0 = 10.\n",
    "concentrations_ = [L0, P0]\n",
    "amplitudes_ = np.array([80, 55, 60])\n",
    "time = np.linspace(0, 20, num=200)\n",
    "noise = 0.02*np.random.randn(len(time))\n",
    "\n",
    "fluorescence_data = forward_model_numpy(rates_, amplitudes_, concentrations_, time) + noise\n",
    "\n",
    "plt.plot(time, fluorescence_data,'o-')\n",
    "plt.plot(time, forward_model_theano(rates_, amplitudes_, concentrations_, time) + noise)\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('fluorescence')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as model:\n",
    "    amplitudes =  pm.Uniform('amplitudes', \n",
    "                             lower=np.array([0., 0., 0.]), \n",
    "                             upper=np.array([100., 100., 100.]), shape=3)\n",
    "    rates =  pm.Uniform('rates', \n",
    "                       np.array([0, 0, 0, 0]), \n",
    "                       np.array([300., 300., 300., 300.]), shape=4)\n",
    "    #rates = pm.Normal('rates', \n",
    "    #                  mu=np.array([1.3, 23.0, 1.3, 6E-4]), \n",
    "    #                  tau=np.array([0.3, 5.0, 0.3, 2E-4])**-2., shape=4)\n",
    "    concentrations = pm.Normal('concentrations', \n",
    "                               mu=np.array([L0, P0]), \n",
    "                               tau=np.array([0.1*L0, 0.1*P0])**-2, shape=2)\n",
    "    \n",
    "    c0 = tt.stacklists([concentrations[1], 0.0, 0.0])\n",
    "    M = tt.stacklists([\n",
    "        [-rates[0]*concentrations[0],  rates[1], 0.0],\n",
    "        [ rates[0]*concentrations[0], -rates[1] - rates[2], rates[3]],\n",
    "        [ 0.0,                         rates[2],           -rates[3]]\n",
    "    ])\n",
    "    eigenvalues, V = tt.nlinalg.eig(M)\n",
    "    V_inv = tt.nlinalg.matrix_inverse(V)\n",
    "    z = tt.dot(V_inv, c0)\n",
    "    a = tt.dot(V.T, amplitudes)\n",
    "    ifmodel = tt.sum(tt.exp(eigenvalues*time[:, np.newaxis])*a*z, axis=1)   \n",
    "\n",
    "    err = pm.Uniform('err', 0, 1000)\n",
    "    mu = pm.Deterministic('mu', ifmodel)\n",
    "    pm.Normal('fluorescence', mu=mu, tau=err, observed=fluorescence_data)\n",
    "    trace = pm.sample(100000, tune=10000) # TODO: increase number of samples!!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot validations\n",
    "stride = 1000\n",
    "for r, a, c in zip(trace['rates'][::stride], \n",
    "                   trace['amplitudes'][::stride], \n",
    "                   trace['concentrations'][::stride]):\n",
    "    plt.plot(time, forward_model_numpy(r, a, c, time), c='cyan', alpha=0.1)\n",
    "plt.plot(time, fluorescence_data, 'ok')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('fluorescence')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(3)\n",
    "axes[0].hist(trace['amplitudes'][:, 0])\n",
    "axes[0].axvline(x=amplitudes_[0], ymin=0, ymax=1000, c='orange')\n",
    "axes[1].hist(trace['amplitudes'][:, 1])\n",
    "axes[1].axvline(x=amplitudes_[1], ymin=0, ymax=1000, c='orange')\n",
    "axes[2].hist(trace['amplitudes'][:, 2])\n",
    "axes[2].axvline(x=amplitudes_[2], ymin=0, ymax=1000, c='orange')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2)\n",
    "axes[0].hist(trace['concentrations'][:, 0])\n",
    "axes[0].axvline(x=concentrations_[0], ymin=0, ymax=1000, c='orange')\n",
    "axes[1].hist(trace['concentrations'][:, 1])\n",
    "axes[1].axvline(x=concentrations_[1], ymin=0, ymax=1000, c='orange')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(4)\n",
    "axes[0].hist(trace['rates'][:, 0])\n",
    "axes[0].axvline(x=rates_[0], ymin=0, ymax=1000, c='orange')\n",
    "axes[1].hist(trace['rates'][:, 1])\n",
    "axes[1].axvline(x=rates_[1], ymin=0, ymax=1000, c='orange')\n",
    "axes[2].hist(trace['rates'][:, 2])\n",
    "axes[2].axvline(x=rates_[2], ymin=0, ymax=1000, c='orange')\n",
    "axes[3].hist(trace['rates'][:, 3])\n",
    "axes[3].axvline(x=rates_[3], ymin=0, ymax=1000, c='orange')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}