{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from lumicks import pylake\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kinesin Attached to a Bead Walking on Microtubule\n",
    "\n",
    "[Download this page as a Jupyter notebook](_downloads/b2c9fbcbc80eebf534af8abb7b9776f3/cytoskeletal_kinesin_bead_open_loop.ipynb)\n",
    "\n",
    "In this assay we had microtubules on the surface. We trapped beads with Kinesin (molecular motor) and had ATP inside the assay. As we lowered the kinesin-coated beads on top of a microtubule, it attached to it and started stepping on them. Kinesins were pulling the bead out of the center of the trap and thus increasing the force on the bead. At one point the kinesins couldn’t keep up with this increased force and the bead construct snapped back to its original position.\n",
    "\n",
    "After this, the cycles starts again with Kinesins pulling the bead out of the center of the trap.\n",
    "\n",
    "With the IRM, you can see unlabeled microtubules and the kinesin-coated bead on top of one of them.\n",
    "\n",
    "To see and characterize this behaviour, we need to have the following plots:\n",
    "\n",
    "* Plot force versus time\n",
    "  \n",
    "* Plot displacement of the bead from the center of the trap over time\n",
    "  \n",
    "Load all the needed libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from lumicks import pylake"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Install Pylake, in case it’s not installed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install lumicks.pylake"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = r'20190215-170635 Marker Kinesin stepping open loop.h5';\n",
    "data = pylake.File(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the contents of the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Force versus Time\n",
    "\n",
    "Get the raw data out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forcey = data['Force HF']['Force 1y']\n",
    "forcex = data['Force HF']['Force 1x']\n",
    "\n",
    "# We need to convert the time from nanoseconds to milliseconds and make sure that time=0 at 0 seconds\n",
    "time = (forcey.timestamps - forcey.timestamps[0]) /1e9\n",
    "\n",
    "forcex_data = forcex.data\n",
    "forcey_data = forcey.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot force in x and y:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(13, 5))\n",
    "\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(time, forcex_data)\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Force X (pN)')\n",
    "\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(time, forcey_data)\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Force Y (pN)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can clearly see that the bead was moving in the y direction, so for now we’re just going to work with that. Later I have an example of how to deal with a bead moving at an angle, like at 45 degrees.\n",
    "\n",
    "But for now, let’s also downsample the force data to 100 Hz and plot the two together.\n",
    "\n",
    "Downsample the y force data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "downsampled_rate = 100 # Hz\n",
    "\n",
    "sample_rate = forcey.sample_rate\n",
    "\n",
    "forcey_downsamp = forcey.downsampled_by(int(sample_rate/downsampled_rate))\n",
    "forcex_downsamp = forcex.downsampled_by(int(sample_rate/downsampled_rate))\n",
    "time_downsampled = (forcey_downsamp.timestamps - forcey_downsamp.timestamps[0]) /1e9\n",
    "\n",
    "forcey_downsamp_data = forcey_downsamp.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two sampling rates are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Original sampling rate is ' + str(sample_rate) + ' Hz')\n",
    "print('Downsampled rate is ' + str(downsampled_rate) + ' Hz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the original force and the downsampled rate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(13, 5))\n",
    "\n",
    "plt.plot(time, forcey_data,label='Original, 30 kHz')\n",
    "plt.plot(time_downsampled, forcey_downsamp_data, 'r',label='Downsampled, 100 Hz')\n",
    "\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Force X (pN)')\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Displacement versus Time\n",
    "\n",
    "We need to convert the force to displacement, which we can do with the following formula:\n",
    "\n",
    "$$Delta x = frac{F}{k}$$\n",
    "\n",
    "where $F$ is the force and $k$ is the trap stiffness. Force we already have, we need to get stiffness.\n",
    "\n",
    "Get stiffness from force calibration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = data['Calibration']['9']['Force 1y'].h5\n",
    "ky = params.attrs.get(\"kappa (pN/nm)\")\n",
    "\n",
    "params = data['Calibration']['9']['Force 1x'].h5\n",
    "kx = params.attrs.get(\"kappa (pN/nm)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The stiffness values are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ky) # this is in pN/nm\n",
    "print(kx) # this is in pN/nm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate and plot displacement versus time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "displacement = forcey_data / ky\n",
    "displacement_downsampled = forcey_downsamp_data / ky\n",
    "\n",
    "\n",
    "plt.figure(figsize=(13, 5))\n",
    "\n",
    "plt.plot(time, displacement,label='Original, 30 kHz')\n",
    "plt.plot(time_downsampled, displacement_downsampled, 'r',label='Downsampled, 100 Hz')\n",
    "\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Displacement (nm)')\n",
    "plt.legend()\n",
    "\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distance and Force versus Time on Same Graph\n",
    "\n",
    "Plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots(figsize=(13, 5))\n",
    "\n",
    "plt.plot(time, displacement,label='Original, 30 kHz')\n",
    "\n",
    "ax1.set_xlabel('Time (s)')\n",
    "ax1.set_ylabel('Displacement (nm)')\n",
    "ax1.set_yticks([-60,-50,-40,-30,-20,-10,0,10,20,30,40,50,60,70,80,90,100])\n",
    "ax1.grid()\n",
    "\n",
    "\n",
    "# create another axis\n",
    "ax2 = ax1.twinx()\n",
    "\n",
    "# ax2.plot(time_downsampled, fy_downsamp.data+5*ky, 'r-')\n",
    "ax2.plot(time_downsampled, forcey_downsamp_data, 'r',label='Downsampled, 100 Hz')\n",
    "\n",
    "ax2.set_ylabel('Force (pN)', color='r')\n",
    "ax2.tick_params('y', colors='r')\n",
    "\n",
    "\n",
    "# Here we just make sure that both the displacement and the force axis have the same limits\n",
    "ylimits = [-60, 100]\n",
    "ylim2 =[]\n",
    "for i in ylimits:\n",
    "    ylim2.append(i*ky)\n",
    "\n",
    "ax1.set_ylim(ylimits)\n",
    "ax2.set_ylim(ylim2)\n",
    "ax1.set_xlim([0, 5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## X vs Y Position of the Bead\n",
    "\n",
    "To get an idea in which direction the microtubule was oriented, which direction the force was applied, we plot the (x,y) position of the bead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(forcex_downsamp.data / kx , forcey_downsamp_data / ky,'.')\n",
    "plt.xlim([-60, 80])\n",
    "plt.ylim([-60, 80])\n",
    "\n",
    "plt.ylabel('y-position (nm)')\n",
    "plt.xlabel('x-position (nm)')\n",
    "plt.grid()"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 2
}