P2X models.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:b612a0ff69e161943730c6e76e22618e0eba7ab857f4469f7454645308b3944b"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {
      "heading_collapsed": false
     },
     "source": [
      "Yan et al. 2010 P2X Markov Model"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {
      "heading_collapsed": false
     },
     "source": [
      "Set variables according to distributions"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.integrate import odeint\n",
      "from numpy import linspace,logspace, asarray, array, concatenate, set_printoptions\n",
      "from itertools import izip_longest\n",
      "from matplotlib import pyplot as plt\n",
      "import itertools\n",
      "import pandas as pd\n",
      "from matplotlib import pylab as pyl\n",
      "import brian"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": [
      {
       "ename": "ImportError",
       "evalue": "No module named brian",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mImportError\u001b[0m                               Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-1-e6b0fec1e6df>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      6\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mpandas\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mpd\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      7\u001b[0m \u001b[1;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mpylab\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mpyl\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 8\u001b[1;33m \u001b[1;32mimport\u001b[0m \u001b[0mbrian\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
        "\u001b[1;31mImportError\u001b[0m: No module named brian"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {
      "heading_collapsed": false
     },
     "source": [
      "Set variables to exact values"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#NOTE: variables below are in units of seconds and Molars!\n",
      "#be aware of this when choosing the value of the ATP concentration!\n",
      "k1 = 0.3\n",
      "k2 = 40000.\n",
      "k3 = 2.4\n",
      "k4 = 50000.\n",
      "k5 = 1.58\n",
      "k6 = 7000\n",
      "\n",
      "L1 = 0.0001\n",
      "L2 = 0.004\n",
      "L3 = 0.5\n",
      "#L3 = 0.1\n",
      "\n",
      "g12 = 1.5e-8\n",
      "g34 = 4.5e-8\n",
      "\n",
      "V = -60.0\n",
      "E = 0.0"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {
      "heading_collapsed": false
     },
     "source": [
      "Piecewise P2X7 Solver"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "set_printoptions(3)\n",
      "\n",
      "def dState(y,t,atp):\n",
      "\n",
      "    q1 = y[0]\n",
      "    q2 = y[1]\n",
      "    q3 = y[2]\n",
      "    q4 = y[3]\n",
      "    c1 = y[4]\n",
      "    c2 = y[5]\n",
      "    c3 = y[6]\n",
      "    c4 = y[7]\n",
      "    \n",
      "    A = atp\n",
      "\n",
      "    dc1 = k1*c2 + L1*c4 - 3.*k2*A*c1\n",
      "    dc2 = 3.*k2*A*c1 + 2.*k3*q1 - (k1+2.*k4*A)*c2\n",
      "    \n",
      "    dq1 = 2.*k4*A*c2 + 3.*k5*q2 - (2.*k3 + k6*A)*q1\n",
      "    dq2 = k6*A*q1 + L2*q3 - (3.*k5 + L3)*q2\n",
      "    dq3 = k2*A*q4 + L3*q2 -(3.*k1 + L2)*q3\n",
      "    dq4 = 2.*k2*A*c3 + 3.*k1*q3 - (2.*k1 + k2*A)*q4\n",
      "    \n",
      "    dc3 = 3*k2*A*c4 + 2.*k1*q4 - (k1 + 2.*k2*A)*c3\n",
      "    dc4 = k1*c3 - (L1 + 3.*k2*A)*c4\n",
      "        \n",
      "    dy = [dq1,dq2,dq3,dq4,dc1,dc2,dc3,dc4]\n",
      "    return dy\n",
      "\n",
      "\n",
      "def single_solve(tStart=0., tEnd=100., dt=.1, doseConcentration=3.2, startConditions=[0., 0., 0., 0., 1., 0., 0., 0.]):\n",
      "    '''\n",
      "    Create time array and return the results of a single run of odeint.\n",
      "    '''\n",
      "    tLen = (tEnd-tStart)/dt\n",
      "    #solve for q vals, then I vals\n",
      "    time = linspace(tStart,tEnd,tLen)\n",
      "    y = odeint(dState, startConditions, time, args=(doseConcentration,))\n",
      "    return y, time\n",
      "    \n",
      "    \n",
      "def multi_solve(durationDosePairs=[(120,3.2)], y0=[0., 0., 0., 0., 1., 0., 0., 0.]): \n",
      "    \n",
      "    stateNames = ['q1','q2','q3','q4','c1','c2','c3','c4']\n",
      "    allStates = []\n",
      "    \n",
      "    #use same dt for each run (for now. May be a good idea to add as an argument later.)\n",
      "    dt = 0.5\n",
      "    \n",
      "    #set starting conditions for first run\n",
      "    startCond = y0\n",
      "    startTime = 0.0\n",
      "    \n",
      "    for pair in durationDosePairs:\n",
      "        \n",
      "        duration, dose = pair\n",
      "        endTime = startTime + duration\n",
      "        \n",
      "        y,t = single_solve(startTime, endTime, dt, dose, startCond)\n",
      "        s = pd.DataFrame(data=y, index=t, columns=stateNames)\n",
      "        #print \"time range: {} - {}:\".format(t[0],t[-1])\n",
      "        \n",
      "        #reset starting time and conditions base on end time and conditions of last run\n",
      "        startCond = y[-1]\n",
      "        startTime = endTime\n",
      "        \n",
      "        allStates.append(s)\n",
      "            \n",
      "    return pd.concat(allStates)\n",
      "\n",
      "\n",
      "def current(q1,q2,q3,q4,Vm=-60):\n",
      "    i12 = g12*(q1+q2)*(Vm-E)\n",
      "    i34 = g34*(q3+q4)*(Vm-E)\n",
      "    I = i12 + i34\n",
      "    return I,i12,i34    \n",
      "    \n",
      "def plotSolution(s,n):\n",
      "\n",
      "    f = pyl.figure(n)\n",
      "        \n",
      "    Iab, Ia ,Ib = current(s['q1'],s['q2'],s['q3'],s['q4'])\n",
      "    sq12 = s['q1']+s['q2']\n",
      "\n",
      "    r,c = 3,1\n",
      "\n",
      "    plt.subplot(r,c,1)\n",
      "    #plt.title(\"Q1 + Q2, C4\")\n",
      "    q12 = sq12.plot(label=\"Q1 + Q2\",legend=True)\n",
      "    c4 = s['c4'].plot(label='C4',legend=True)\n",
      "\n",
      "    plt.subplot(r,c,2)\n",
      "    #plt.title(\"Q3, Q4\")\n",
      "    s['q3'].plot(label=\"Q3\",legend=True)\n",
      "    s['q4'].plot(label=\"Q4\",legend=True)\n",
      "\n",
      "    plt.subplot(r,c,3)\n",
      "    #plt.title(\"I1, I2\")\n",
      "    Ia.plot(label=\"I1\",legend=True)\n",
      "    Ib.plot(label=\"I2\",legend=True)\n",
      "    Iab.plot(label=\"I1+I2\",legend=True)\n",
      "    plt.tight_layout()\n",
      "    \n",
      "    return f\n",
      "    #plt.show()    "
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#ATP values (but const params still derived from BzATP data)\n",
      "m1 = (100*10**-3)*55 #Molars\n",
      "s1 = multi_solve([(10.,0.0),(120.,m1),(120.,0.0)])\n",
      "#BzATP\n",
      "m2= (100*10**-6)#(3.2*10**-6)\n",
      "s2 = multi_solve([(10.,0.0),(120.,m2),(120.,0.0)])\n",
      "\n",
      "Iab= current(s1['q1'],s1['q2'],s1['q3'],s1['q4'])[0]\n",
      "Iab2= current(s2['q1'],s2['q2'],s2['q3'],s2['q4'])[0]\n",
      "\n",
      "'''\n",
      "f = pyl.figure(0)\n",
      "plt.subplot(1,2,1)\n",
      "Iab.plot(label=\"I1+I2\",legend=True)\n",
      "\n",
      "plt.subplot(1,2,2)\n",
      "Iab2.plot(label=\"I1+I2\",legend=True)\n",
      "'''\n",
      "\n",
      "f, (ax1, ax2) = plt.subplots(1, 2, sharey=True)\n",
      "ax1.plot(Iab.keys(), Iab.values)\n",
      "ax1.set_title(\"ATP=\"+str(m1))\n",
      "ax2.plot(Iab2.keys(), Iab2.values)\n",
      "ax2.set_title(\"BzATP=\"+str(m2))\n",
      "\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#If you want to set variables to random values, use this rather than the variable definitions up top.\n",
      "\n",
      "from numpy import random\n",
      "normDist = random.normal\n",
      "uniDist = random.uniform\n",
      "\n",
      "atp = 32\n",
      "\n",
      "k1 = normDist(0.3,0.03)\n",
      "k2 = normDist(40000.,4000.)\n",
      "k3 = normDist(2.4,.24)\n",
      "k4 = normDist(50000., 5000.)\n",
      "k5 = normDist(1.58,.158)\n",
      "k6 = normDist(7000,700.)\n",
      "\n",
      "L1 = 0.0001\n",
      "L2 = normDist(0.004,4*(10**-4))\n",
      "L3 = uniDist(0.1,3.)\n",
      "\n",
      "g12 = normDist(1.5e10-8,1.5e-9)\n",
      "g34 = normDist(4.5e10-8, 4.5e-9)\n",
      "\n",
      "E = normDist(0.0,1)\n",
      "V = -60.0\n"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "t = np.vstack(ts)\n",
      "#y = np.vstack(ys)\n",
      "print len(ts)\n",
      "#y = odeint(dfv, startConditions, time, args=(doseConcentration,))\n",
      "s = pd.DataFrame(data=ys, index=t, columns=stateNames)\n",
      "s['v'].plot()\n",
      "if len(ts) > 1: plt.show()"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "s1=s.fillna(0)\n",
      "print s1\n",
      "s1['v'].plot()\n",
      "plt.show()\n",
      "\n",
      "#iatp = current(s['q1'],s['q2'],s['q3'],s['q4'],s['v'])[0]\n",
      "iatp = current(s1['q1'],s1['q2'],s1['q3'],s1['q4'])[0]\n",
      "iatp.plot()\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "s=s1\n",
      "\n",
      "Iab, Ia ,Ib = current(s['q1'],s['q2'],s['q3'],s['q4'])\n",
      "Iab.plot(label=\"I1+I2\",legend=True)\n",
      "plt.show()\n",
      "vm = linspace(V,-5,len(s['q1']))\n",
      "Iab, Ia ,Ib = current(s['q1'],s['q2'],s['q3'],s['q4'],vm)\n",
      "Iab.plot(label=\"I1+I2\",legend=True)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {
      "heading_collapsed": false
     },
     "source": [
      "P2X7 Model with Burster"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.integrate as spi\n",
      "import numpy as np\n",
      "#from scipy.integrate import ode\n",
      "import matplotlib.pyplot as plt\n",
      "import warnings\n",
      "import pandas as pd\n",
      "import math\n",
      "#### testing voltage dependence ####\n",
      "\n",
      "#------------------------------------------------------------------------------------------\n",
      "#### INITIAL CONDITIONS ####\n",
      "Vs0=-60.0 #initial cell memberane voltage\n",
      "ns0=0.004 #initial value of gating var n\n",
      "hs0=0.33 #initial value of gating var h\n",
      "Vd0=-50.0 #initial dendritic membrane voltage\n",
      "C0=0.03 #initial calcium 2+ balance\n",
      "l0=0.93 # initial value of IP3 channel gating variable\n",
      "\n",
      "#### TIME (ms) ####\n",
      "# Duration\n",
      "seconds = 10 #set to the number of seconds you want to see\n",
      "# time step\n",
      "dt=0.3 #(ms)\n",
      "\n",
      "#### MEMBRANE CAPACITANCE (pF) ####\n",
      "Cms=21\n",
      "#dendrite membrane capacitance\n",
      "Cmd=5\n",
      "\n",
      "#### REVERSAL POTENTIALS (mV) ####\n",
      "# sodium potential (for normal sodium current and persistent sodium current)\n",
      "ena=50\n",
      "# potassium potential\n",
      "ek=-85\n",
      "# leak reversal potential\n",
      "eL=-54.0 # suggested value for vleak from Butera paper are -60,-57.5,-54\n",
      "# reversal potential of non-NMDA glutamatergic synaptic currents.\n",
      "eSyn=0.0\n",
      "\n",
      "#### CONDUCTANCE (nS) ####\n",
      "# potassium current conductance\n",
      "gk=11.2\n",
      "# sodium current conductance\n",
      "gna=28\n",
      "# persistent sodium current conductance\n",
      "gnaps=2.8\n",
      "# leak current conductance\n",
      "gL=5.6\n",
      "# tonic conductance\n",
      "gtonic=0.0\n",
      "\n",
      "gcan = 1.0\n",
      "\n",
      "#### APPLIED CURRENT ####\n",
      "Iaps=0\n",
      "\n",
      "#### HALF (IN)ACTIVATION VOLTAGE (mV) ####\n",
      "vm=-34\n",
      "vn=-29\n",
      "vmp=-40\n",
      "vh=-48\n",
      "\n",
      "#### INVERSE SLOPE MULTIPLIER (mV) ####\n",
      "sm=-5\n",
      "sn=-4\n",
      "smp=-6\n",
      "sh=6\n",
      "\n",
      "#### TIME CONSTANT MAX VALUES (ms) #####\n",
      "taunb=10 #tau max for n\n",
      "tauhb=10000 #tau max for h\n",
      "\n",
      "\n",
      "#### Constants for calculating Ca Flux: ER --> Cytosol ####\n",
      "I=1 #IP3 concentration\n",
      "LL=0.37 #ER leak permeability\n",
      "P=31000 #maximum total permeability of IP3 channels\n",
      "Ki=1.0 #dissociation consts for IP3 receptor activity by IP3\n",
      "Ka=0.4 #dissociation consts for IP3 receptor activity by Ca\n",
      "\n",
      "#### Constants for calculating Ca Flux: Cytosol --> ER ####\n",
      "Ve=400 #Maximal SERCA pump rate\n",
      "Ke=0.2 #coefficient for SERCA pumps\n",
      "\n",
      "#### ER Ca CONCENTRATION ####\n",
      "Ct=1.25 #Total Ca\n",
      "sigma=0.185 #ratio of cytosolic to ER volume\n",
      "\n",
      "#The ER parameters\n",
      "fi=0.0001 #bound Ca concentration in cytosol\n",
      "Vi=4 #free Ca concentration in cytosol\n",
      "A=0.005 #scaling const.\n",
      "Kd=0.4 #dissociation constant for IP3 receptor inactivation by Ca\n",
      "\n",
      "#### Calcium current activation ####\n",
      "Kcan=0.74 # microM \n",
      "ncan=0.97\n",
      "\n",
      "#### Ratio of somatic to total area ####\n",
      "k=0.3\n",
      "\n",
      "#---------------------------\n",
      "#NOTE: variables below are in units of ms and Molars!\n",
      "#be aware of this when choosing the value of the ATP concentration!\n",
      "k1 = 0.3e-3 \n",
      "k2 = 40000e-3\n",
      "k3 = 2.4e-3\n",
      "k4 = 50000.0e-3\n",
      "k5 = 1.58e-3\n",
      "k6 = 7000e-3\n",
      "L1 = 0.0001e-3\n",
      "L2 = 0.004e-3\n",
      "L3 = 0.5e-3\n",
      "#L3 = 0.1e-3\n",
      "\n",
      "g12 = 15 # nS\n",
      "g34 = 45 # nS\n",
      "E = 0.0 # mV\n",
      "\n",
      "\n",
      "def dfv(y,t,atp):\n",
      "    \n",
      "    q1 = y[0]\n",
      "    q2 = y[1]\n",
      "    q3 = y[2]\n",
      "    q4 = y[3]\n",
      "    c1 = y[4]\n",
      "    c2 = y[5]\n",
      "    c3 = y[6]\n",
      "    c4 = y[7]\n",
      "    Vs = y[8]\n",
      "    ns = y[9]\n",
      "    hs = y[10]\n",
      "    Vd = y[11] #Dendritic membrane voltage\n",
      "    C = y[12] #Intracellular Ca2+ balance\n",
      "    l = y[13] #IP3 channel gating variable\n",
      "    \n",
      "    # SOMATIC FUNCTIONS\n",
      "    minfs = 1/(1+math.exp((Vs-vm) /sm))\n",
      "    ninfs = 1/(1+math.exp((Vs-vn) /sn))\n",
      "    minfps = 1/(1+math.exp((Vs-vmp)/smp))\n",
      "    hinfs = 1/(1+math.exp((Vs-vh) /sh))\n",
      "    \n",
      "    tauns = taunb/math.cosh((Vs-vn)/(2*sn))\n",
      "    tauhs = tauhb/math.cosh((Vs-vh)/(2*sh))\n",
      "    \n",
      "    # DENDRITIC FUNCTIONS\n",
      "    #Calculate ER Ca\n",
      "    Ce = (Ct - C)/sigma\n",
      "    # Flux of Ca from ER to cytosol(regulated by IP3 receptors)\n",
      "    J_ER_in=(LL + P*math.pow(((I*C*l)/((I+Ki)*(C+Ka))),3))*(Ce - C) \n",
      "    # Flux from cytosol back to ER (controlled by SERCA pumps)\n",
      "    J_ER_out=Ve*math.pow(C,2)/(math.pow(Ke,2)+math.pow(C,2))\n",
      "    # Activation of calcium current (I_can)\n",
      "    caninf = 1/(1+math.pow((Kcan/C),ncan))\n",
      "    \n",
      "    I_nas = gna*math.pow(minfs,3)*(1-ns)*(Vs-ena) #sodium current\n",
      "    I_ks = gk*math.pow(ns,4)*(Vs-ek) #potassium current\n",
      "    I_naps = gnaps*minfps*hs*(Vs-ena) #persistent sodium current\n",
      "    I_L = gL*(Vs-eL) #pA\n",
      "    I_atp = g12*(q1+q2)*(1/10)*(Vs) + g34*(q3+q4)*(Vs)\n",
      "    I_ds = gcan*(Vd-Vs)/k # modification of somatic current due to dendritic current\n",
      "    \n",
      "    I_can = gcan*caninf*(Vd-ena) # calcium current\n",
      "    I_sd = gcan*(Vs-Vd)/(1-k) # modification of dendritic current due to somatic current\n",
      "    \n",
      "    #dv = (-iatp-ileak)/C # i is in units of picoamps and its around 10**12 pA, so dv is in units of Volts!\n",
      "    dVs = (-I_ks-I_nas-I_naps-I_L-I_sd-I_atp)/Cms\n",
      "    dns = (ninfs-ns)/tauns\n",
      "    dhs = (hinfs-hs)/tauhs\n",
      "    \n",
      "    # DENDRITIC EQUATIONS\n",
      "    dVd = (-I_can-I_ds)/Cmd\n",
      "    dC = fi/Vi*( J_ER_in - J_ER_out)\n",
      "    dl = A*( Kd - (C + Kd)*l )\n",
      "    \n",
      "    #P2X7 EQUATIONS\n",
      "    dc1 = k1*c2 + L1*c4 - 3.*k2*atp*c1\n",
      "    dc2 = 3.*k2*atp*c1 + 2.*k3*q1 - (k1+2.*k4*atp)*c2\n",
      "    \n",
      "    dq1 = 2.*k4*atp*c2 + 3.*k5*q2 - (2.*k3 + k6*atp)*q1\n",
      "    dq2 = k6*atp*q1 + L2*q3 - (3.*k5 + L3)*q2\n",
      "    dq3 = k2*atp*q4 + L3*q2 -(3.*k1 + L2)*q3\n",
      "    dq4 = 2.*k2*atp*c3 + 3.*k1*q3 - (2.*k1 + k2*atp)*q4\n",
      "    \n",
      "    dc3 = 3*k2*atp*c4 + 2.*k1*q4 - (k1 + 2.*k2*atp)*c3\n",
      "    dc4 = k1*c3 - (L1 + 3.*k2*atp)*c4\n",
      "     \n",
      "    dy = [dq1,dq2,dq3,dq4,dc1,dc2,dc3,dc4,dVs,dns,dhs,dVd,dC,dl]     \n",
      "    return dy\n",
      "\n",
      "\n",
      "stateNames = ['q1','q2','q3','q4','c1','c2','c3','c4','v']\n",
      "\n",
      "tStart=0.\n",
      "tEnd=5000.\n",
      "dt=1.\n",
      "atp=0.#3.2e-6 #Molars #equal to 3.2 mM\n",
      "startConds=[0., 0., 0., 0., 1., 0., 0., 0.,Vs0, ns0, hs0,Vd0,C0,l0]\n",
      "tLen = (tEnd-tStart)/dt\n",
      "time = np.linspace(tStart,tEnd,tLen)\n",
      "try:\n",
      "    y=spi.odeint(dfv,y0=startConds,t=time, args=(atp,))\n",
      "except Exception as e:\n",
      "    print e\n",
      "s = pd.DataFrame(data=y[:,0:9], index=time, columns=stateNames)\n",
      "\n",
      "s['v'].plot()\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The following is an experiment to see what happens when I plot the functions used in Yan et al. 2010 to fit their data."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "from numpy.linalg import norm\n",
      "def phase1(x,tau,A=3.2):\n",
      "    #A=3.2\n",
      "    #tau = 17.\n",
      "    C = 0\n",
      "    return A*(1-math.exp(-x/tau))+C\n",
      "\n",
      "def deact(x):\n",
      "    A1 = 3.2\n",
      "    A2 = 3.2\n",
      "    tau1 = 3\n",
      "    tau2 = 2\n",
      "    C = 3\n",
      "    return A1*math.exp(-x/tau1) + A2*math.exp(-x/tau2) + C\n",
      "\n",
      "def phase2(x):\n",
      "    return\n",
      "def plot(f):\n",
      "    tStart = 0.0\n",
      "    tEnd=50.\n",
      "    dt=.1\n",
      "    t = linspace(tStart,tEnd,tEnd/dt)\n",
      "    plt.plot(t, [f(x) for x in t])\n",
      "    plt.show()\n",
      "    \n",
      "#plot(phase1, 64.8)\n",
      "\n",
      "f=phase1\n",
      "tStart = 0.0\n",
      "tEnd=50.\n",
      "dt=.1\n",
      "norm =100\n",
      "t = linspace(tStart,tEnd,tEnd/dt)\n",
      "tau_1=64.8*1000. #ms\n",
      "tau_3=78.6*1000. #ms\n",
      "tau_4=5.6 #seconds\n",
      "tau_2=10. #seconds\n",
      "c_13 = 10\n",
      "c_24=100\n",
      "p2x1 = asarray([f(x,tau_1,c_13) for x in t])\n",
      "p2x3 = asarray([f(x,tau_3,c_13) for x in t])\n",
      "p2x2 = asarray([f(x,tau_2,c_24) for x in t])\n",
      "p2x4 = asarray([f(x,tau_4,c_24) for x in t])\n",
      "\n",
      "#plt.plot(t, p2x1)\n",
      "#plt.plot(t, (p2x3))\n",
      "#plt.plot(t, (p2x2))\n",
      "plt.plot(t, (p2x4))\n",
      "plt.plot()\n",
      "plt.show()\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "#plot(deact)"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy import *\n",
      "import matplotlib.pyplot as plt\n",
      "start = -.5\n",
      "end = 4\n",
      "atpRange = logspace(start,end)\n",
      "ec50 = 14.5\n",
      "h = 0.9\n",
      "\n",
      "def i_atp(atp):\n",
      "    return (1 + (ec50/atp)**h)**-1\n",
      "\n",
      "plt.semilogx(atpRange, [i_atp(atp) for atp in atpRange])\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Riedel et al 2007\n",
      "\n",
      "g = 9.7e-12\n",
      "E = -(1.3e-3)\n",
      "V = linspace(-100,60)\n",
      "\n",
      "plt.plot(V, g*V - E*g)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Riedel model\n",
      "alpha = 200.\n",
      "beta = 66.\n",
      "k1 = 50.\n",
      "k2 = 360.\n",
      "atp = 3.2*10\n",
      "\n",
      "Vm=-60\n",
      "g=9.5\n",
      "\n",
      "def riedelModel(y,t,args=None):\n",
      "    \n",
      "    A = y[0]\n",
      "    B = y[1]\n",
      "    C = y[2]\n",
      "    E = y[3]\n",
      "    \n",
      "    dA = k2*B - 2*k1*atp*A\n",
      "    dB = -(k1*atp + k2)*B + 2*k1*atp*A + 2*k2*C\n",
      "    dC = -(beta + 2*k2)*C + alpha*E +k1*atp*B\n",
      "    dE = beta*C - alpha*E\n",
      "    \n",
      "    return [dA,dB,dC,dE]\n",
      "\n",
      "\n",
      "y0 = asarray([1,0,0,0])\n",
      "s = 0.\n",
      "e = 200.\n",
      "dt = .5\n",
      "t = linspace(s,e,(e-s)/dt)\n",
      "y1 = odeint(riedelModel,y0,t,args=(atp,))\n",
      "y=y1.T\n",
      "plt.plot(t,y[0])  \n",
      "plt.plot(t,y[1]) \n",
      "plt.show()\n",
      "plt.plot(t,y[2])  \n",
      "plt.plot(t,y[3])  \n",
      "plt.show()\n",
      "I=[g * (Vm-0) * O for O in y1[:,3]]\n",
      "plt.plot(t,I)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "#Reidel\n",
      "kd = .117\n",
      "n = 1.4\n",
      "P0inf = 0.26\n",
      "t_start = .01\n",
      "t_end = 10**3\n",
      "dt = .1\n",
      "atp = logspace(-2,3)\n",
      "P0 = map(lambda atp: P0inf/(1+(10**-math.log10(kd))/atp)**n, atp)\n",
      "plt.xscale('log')\n",
      "\n",
      "reidelplot = plt.plot(atp, P0)\n",
      "#reidelplot\n",
      "plt.show()\n",
      "#Bean\n",
      "\n",
      "i_bean = map(lambda atp : 1.04/(1+(2.7/atp)),atp)\n",
      "beanplot = plt.plot(atp, i_bean)\n",
      "plt.show()\n",
      "\n"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from urlparse import urlparse\n",
      "import sys\n",
      "\n",
      "#import gtk\n",
      "#import tkSimpleDialog\n",
      "import matplotlib\n",
      "import matplotlib.image as mpimg\n",
      "import matplotlib.pyplot as plt\n",
      " \n",
      "import numpy as np\n",
      "    \n",
      "def tellme(s):\n",
      "    print s\n",
      "    sys.stdout.flush()\n",
      "    plt.title(s,fontsize=16)\n",
      "    plt.draw()\n",
      " \n",
      "def pic2data(source,straight=True):\n",
      "    \"\"\" GUI to get data from a XY graph image. Either provide the graph\n",
      "        as a path to an image in 'source' or copy it to the clipboard.\n",
      "    \"\"\"\n",
      "     \n",
      "    ##### GET THE IMAGE\n",
      "    image = mpimg.imread(source)\n",
      "    origin = 'lower'\n",
      " \n",
      "    ###### DISPLAY THE IMAGE\n",
      "\n",
      "    plt.ion() # interactive mode !\n",
      "    fig, ax = plt.subplots(1)\n",
      "    imgplot = ax.imshow(image, origin=origin)\n",
      "    fig.canvas.draw()\n",
      "    plt.draw()\n",
      "     \n",
      "    ##### PROMPT THE AXES\n",
      "     \n",
      "    def promptPoint(text=None):\n",
      "         \n",
      "        if text is not None: tellme(text)\n",
      "        return  np.array(plt.ginput(1,timeout=-1)[0])\n",
      "     \n",
      "    def askValue(text='',initialvalue=0.0):\n",
      "        a = float(raw_input(text))\n",
      "        return a\n",
      "        '''\n",
      "        return tkSimpleDialog.askfloat(text, 'Value:',\n",
      "                                         initialvalue=initialvalue)\n",
      "        '''\n",
      "    origin = promptPoint('Place the origin ')\n",
      "    origin_value = askValue('X origin ',0),askValue('Y origin ',0)\n",
      "                   \n",
      "    Xref =  promptPoint('Place the X reference ')\n",
      "    Xref_value = askValue('X reference ',1.0)\n",
      "     \n",
      "    Yref =  promptPoint('Place the Y reference ')\n",
      "    Yref_value = askValue('Y reference ',1.0)\n",
      "    if straight :\n",
      "         \n",
      "        Xref[1] = origin[1]\n",
      "        Yref[0] = origin[0]\n",
      "\n",
      "    ##### PROMPT THE POINTS \n",
      "    selected_points = []\n",
      "     \n",
      "    tellme(\"Select your points !\")\n",
      "    print \"Right-click or press 's' to select\"\n",
      "    print \"Left-click or press 'del' to deselect\"\n",
      "    print \"Middle-click or press 'Enter' to confirm\"\n",
      "    print \"Note that the keyboard may not work.\"\n",
      "    sys.stdout.flush()\n",
      "\n",
      "    selected_points = plt.ginput(-1,timeout=-1)\n",
      "    print \"points: \",selected_points\n",
      "    ##### RETURN THE POINTS COORDINATES\n",
      "     \n",
      "    #~ selected_points.sort() # sorts the points in increasing x order\n",
      "     \n",
      "    # compute the coordinates of the points in the user-defined system\n",
      "     \n",
      "    OXref = Xref - origin\n",
      "    OYref = Yref - origin\n",
      "    xScale =  (Xref_value - origin_value[0]) / np.linalg.norm(OXref)\n",
      "    print \"OXref: \",OXref\n",
      "    print \"xScale: \",xScale\n",
      "#    xScale =  (Xref_value - origin_value[0]) / np.linalg.norm(OXref)\n",
      "\n",
      "    yScale =  (Yref_value - origin_value[1]) / np.linalg.norm(OYref)\n",
      "    \n",
      "    ux = OXref / np.linalg.norm(OXref)\n",
      "    print \"ux: \", ux\n",
      "    uy = OYref / np.linalg.norm(OYref)\n",
      "    \n",
      "    result = [(ux.dot(pt - origin) * xScale + origin_value[0],\n",
      "               uy.dot(pt - origin) * yScale + origin_value[1])\n",
      "               for pt in selected_points ]\n",
      "     \n",
      "    plt.ioff()\n",
      "    print locals()\n",
      "    \n",
      "    return result\n",
      "\n"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "points = pic2data('/Users/morganfine-morris/Desktop/BzATP vs ATP dose response.png')"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Place the origin \n"
       ]
      },
      {
       "name": "stdout",
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "X origin .001\n"
       ]
      },
      {
       "name": "stdout",
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Y origin 0\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Place the X reference \n"
       ]
      },
      {
       "name": "stdout",
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "X reference 10\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Place the Y reference \n"
       ]
      },
      {
       "name": "stdout",
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Y reference 5\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Select your points !\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Right-click or press 's' to select\n",
        "Left-click or press 'del' to deselect\n",
        "Middle-click or press 'Enter' to confirm\n",
        "Note that the keyboard may not work.\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "points:  [(96.612903225806434, 32.383064516129025), (144.27016129032256, 34.181451612903217), (193.72580645161287, 71.947580645161281), (240.48387096774195, 151.9758064516129), (288.14112903225805, 178.05241935483872), (240.48387096774195, 34.181451612903217), (288.14112903225805, 35.979838709677409), (336.69758064516128, 44.072580645161281), (384.35483870967738, 139.38709677419354), (432.91129032258061, 212.22177419354836)]\n",
        "OXref:  [ 386.65322581    0.        ]"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "xScale:  0.0258603816874\n",
        "ux:  [ 1.  0.]\n",
        "{'origin': array([ 48.05645161,  22.49193548]), 'Xref': array([ 434.70967742,   22.49193548]), 'image': array([[[ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        ..., \n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.]],\n",
        "\n",
        "       [[ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        ..., \n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.]],\n",
        "\n",
        "       [[ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        ..., \n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.]],\n",
        "\n",
        "       ..., \n",
        "       [[ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        ..., \n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.]],\n",
        "\n",
        "       [[ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        ..., \n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.]],\n",
        "\n",
        "       [[ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        ..., \n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.],\n",
        "        [ 1.,  1.,  1.,  1.]]], dtype=float32), 'promptPoint': <function promptPoint at 0x10ac00a28>, 'result': [(1.2566883720930231, 0.24886877828054302), (2.4891232558139533, 0.29411764705882354), (3.768065116279069, 1.244343891402715), (4.9772465116279081, 3.2579185520362), (6.2096813953488379, 3.9140271493212677), (4.9772465116279081, 0.29411764705882354), (6.2096813953488379, 0.33936651583710403), (7.4653697674418611, 0.54298642533936659), (8.6978046511627909, 2.9411764705882355), (9.9534930232558132, 4.7737556561085972)], 'origin_value': (0.001, 0.0), 'ax': <matplotlib.axes._subplots.AxesSubplot object at 0x102472790>, 'pt': (432.91129032258061, 212.22177419354836), 'yScale': 0.025160805957429543, 'Xref_value': 10.0, 'source': '/Users/morganfine-morris/Desktop/BzATP vs ATP dose response.png', 'fig': <matplotlib.figure.Figure object at 0x1024670d0>, 'selected_points': [(96.612903225806434, 32.383064516129025), (144.27016129032256, 34.181451612903217), (193.72580645161287, 71.947580645161281), (240.48387096774195, 151.9758064516129), (288.14112903225805, 178.05241935483872), (240.48387096774195, 34.181451612903217), (288.14112903225805, 35.979838709677409), (336.69758064516128, 44.072580645161281), (384.35483870967738, 139.38709677419354), (432.91129032258061, 212.22177419354836)], 'Yref_value': 5.0, 'askValue': <function askValue at 0x10ac006e0>, 'OXref': array([ 386.65322581,    0.        ]), 'imgplot': <matplotlib.image.AxesImage object at 0x10ac90b50>, 'Yref': array([  48.05645161,  221.21370968]), 'xScale': 0.025860381687350092, 'uy': array([ 0.,  1.]), 'ux': array([ 1.,  0.]), 'straight': True, 'OYref': array([   0.        ,  198.72177419])}\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pts=np.array(points)\n",
      "x=pts[:,0]\n",
      "y=pts[:,1]\n",
      "print x\n",
      "print \n",
      "print y\n"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 1.25668837  2.48912326  3.76806512  4.97724651  6.2096814   4.97724651\n",
        "  6.2096814   7.46536977  8.69780465  9.95349302]\n",
        "\n",
        "[ 0.24886878  0.29411765  1.24434389  3.25791855  3.91402715  0.29411765\n",
        "  0.33936652  0.54298643  2.94117647  4.77375566]\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "'origin': array([ 48.05645161,  23.39112903]), \n",
      "'Xref': array([ 435.60887097,   23.39112903]), \n",
      "'image': array([[[ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        ..., \n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.]],\n",
      "\n",
      "       [[ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        ..., \n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.]],\n",
      "\n",
      "       [[ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        ..., \n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.]],\n",
      "\n",
      "       ..., \n",
      "       [[ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        ..., \n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.]],\n",
      "\n",
      "       [[ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        ..., \n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.]],\n",
      "\n",
      "       [[ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        ..., \n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.],\n",
      "        [ 1.,  1.,  1.,  1.]]], dtype=float32), \n",
      "        'result': [(1.2305754060324829, 0.24886877828054293), (2.5065498839907194, 0.29411764705882382), (3.7361252900232018, 1.1990950226244341), (4.965700696055686, 3.1900452488687776), (6.2184756380510446, 3.8914027149321266)], \n",
      "        'origin_value': (0.001, 0.0), \n",
      "        'ax': <matplotlib.axes._subplots.AxesSubplot object at 0x111b77c10>, \n",
      "        'pt': (289.04032258064512, 178.05241935483872), \n",
      "        'yScale': 0.025160805957429536, \n",
      "        'Xref_value': 10.0, \n",
      "        'selected_points': [(95.713709677419359, 33.282258064516114), (145.16935483870967, 35.08064516129032), (192.82661290322579, 71.048387096774178), (240.48387096774195, 150.17741935483869), (289.04032258064512, 178.05241935483872)], \n",
      "        'Yref_value': 5.0, \n",
      "        'OXref': array([ 387.55241935,    0.        ]), \n",
      "        'Yref': array([  48.05645161,  222.11290323]), \n",
      "        'xScale': 0.025800380801764593, \n",
      "        'uy': array([ 0.,  1.]), \n",
      "        'ux': array([ 1.,  0.]), \n",
      "        'OYref': array([   0.        ,  198.72177419])}"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#?curve_fit\n",
      "#fitting the sigmoid functions from fig 1 of Yan et al\n",
      "from scipy.optimize import curve_fit\n",
      "from numpy import polyfit, set_printoptions, logspace, array\n",
      "\n",
      "points = [(96.612903225806434, 32.383064516129025), (142.47177419354838, 35.08064516129032), \\\n",
      "          (192.82661290322579, 71.947580645161281), (239.58467741935482, 151.9758064516129), \\\n",
      "          (288.14112903225805, 178.95161290322579), (240.48387096774195, 33.282258064516114), \\\n",
      "          (288.14112903225805, 35.979838709677409), (337.59677419354836, 44.072580645161281), \\\n",
      "          (385.25403225806451, 138.48790322580643), (433.81048387096774, 211.32258064516128)]\n",
      "\n",
      "px,py = zip(*points)\n",
      "plt.plot(px,py,'o')\n",
      "plt.show()\n",
      "\n",
      "def sigmoid(atp,EC50, h, norm):\n",
      "    return norm*(1+(EC50/atp)**h)**-1\n",
      "\n",
      "#data to fit to\n",
      "bz = (array([.0032,.01,.032,.1,.32]),array([0.4167,0.417,1.2,3.125,3.75]))\n",
      "atp = ([.1,.32,1.,3.2,10.],[0.4167,0.417,0.625, 3.,4.75])\n",
      "\n",
      "#fit data\n",
      "bz_params, other = curve_fit(sigmoid,bz[0],bz[1],[5.1*10**-2,1,1])\n",
      "atp_params, other = curve_fit(sigmoid,atp[0],atp[1],[2.88,1,1])\n",
      "#print param vals\n",
      "print \"BzATP: EC50 = {}, h = {}, norm = {}\".format(*bz_params)\n",
      "print \"ATP: EC50 = {}, h = {}, norm = {}\".format(*atp_params)\n",
      "\n",
      "x=logspace(-3,1.5)\n",
      "plt.xscale('log')\n",
      "\n",
      "#plot points\n",
      "plt.plot(bz[0],bz[1],'o',color=\"k\")\n",
      "plt.plot(atp[0],atp[1],'o',color='r')\n",
      "\n",
      "plt.plot(x,sigmoid(x, *bz_params))\n",
      "plt.plot(x,sigmoid(x, *atp_params))\n",
      "\n",
      "plt.xlim(xmax = 14.)\n",
      "\n",
      "plt.show()\n",
      "\n"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "BzATP: EC50 = 0.0499577838057, h = 1.53734542805, norm = 4.02968947368\n",
        "ATP: EC50 = 2.80137758242, h = 1.71533737463, norm = 5.30182428632\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#not this one\n",
      "import numpy as np\n",
      "from matplotlib import pyplot as plt\n",
      "from matplotlib import pylab\n",
      "from scipy.optimize import curve_fit\n",
      "points = [(96.612903225806434, 32.383064516129025), (142.47177419354838, 35.08064516129032), \\\n",
      "          (192.82661290322579, 71.947580645161281), (239.58467741935482, 151.9758064516129), \\\n",
      "          (288.14112903225805, 178.95161290322579), (240.48387096774195, 33.282258064516114), \\\n",
      "          (288.14112903225805, 35.979838709677409), (337.59677419354836, 44.072580645161281), \\\n",
      "          (385.25403225806451, 138.48790322580643), (433.81048387096774, 211.32258064516128)]\n",
      "\n",
      "px,py = zip(*points)\n",
      "plt.plot(px,py,'o')\n",
      "plt.show()\n",
      "    \n",
      "    \n",
      "def sigmoid(x, x0, k):\n",
      "     y = 1 / (1 + np.exp(-k*(x-x0)))\n",
      "     return y\n",
      " \n",
      "bzatp = ([.0032,.01,.032,.1,.32],[0.4167,0.417,1.2,3.125,3.75])\n",
      "atp = ([.1,.32,1.,3.2,10.],[0.4167,0.417,0.625, 3.,4.75])\n",
      " \n",
      "popt, pcov = curve_fit(sigmoid, bzatp[0], bzatp[1])\n",
      "print popt\n",
      " \n",
      "x_array = np.linspace(.001, 10, 50)\n",
      " = sigmoid(x, *popt)\n",
      " \n",
      "plt.xscale('log')\n",
      "pylab.plot(bzatp[0], bzatp[1], 'o', label='data')\n",
      "pylab.plot(x_array,y, label='fit')\n",
      "pylab.ylim(0, 1.05)\n",
      "pylab.legend(loc='best')\n",
      "pylab.show()"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": [
      {
       "ename": "IndentationError",
       "evalue": "unexpected indent (<ipython-input-3-05ac24d8ed41>, line 24)",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-3-05ac24d8ed41>\"\u001b[1;36m, line \u001b[1;32m24\u001b[0m\n\u001b[1;33m    = sigmoid(x, *popt)\u001b[0m\n\u001b[1;37m    ^\u001b[0m\n\u001b[1;31mIndentationError\u001b[0m\u001b[1;31m:\u001b[0m unexpected indent\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {
      "heading_collapsed": false
     },
     "source": [
      "Curve Fitting Dose-Response Data for Logged and Unlogged X Values"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "proof of concept/how-to"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#the following shows both methods\n",
      "import numpy as np\n",
      "from matplotlib import pyplot as plt\n",
      "from scipy.optimize import curve_fit\n",
      "from numpy import polyfit, set_printoptions\n",
      "import sys\n",
      "\n",
      "#two different fitting functions that produce the same result\n",
      "def Log10TransXVals(x,a1,a2,logx10, p):\n",
      "    '''\n",
      "    y = (a1+(a2-a1)/(1+10**((logx10-x)*p)))\n",
      "    '''\n",
      "    y = (a1+(a2-a1)/(1+10**((logx10-x)*p)))\n",
      "    return y\n",
      "\n",
      "def NonLogXVals(x, a1, a2, ec50, h):\n",
      "    '''\n",
      "    y = a2 +((a1-a2)/(1+(x/ec50)**-h))\n",
      "    '''\n",
      "    y = a2 +((a1-a2)/(1+(x/ec50)**-h))\n",
      "    return y\n",
      "\n",
      "#test datapoints\n",
      "xarr = np.array([-8,-7,-6.5,-6,-5.5,-5,-4]) #exponents, for use with Log10TransXVals\n",
      "xarr_10 = 10**(xarr) #actual values, for use with NonLogXVals\n",
      "yarr = [0,5,20,45,75,95,100] # y coords for use with both functions\n",
      "\n",
      "\n",
      "############################\n",
      "#using Log10TransXVals equation for fitting\n",
      "\n",
      "po = curve_fit(Log10TransXVals, xarr, yarr, [.49, 100.62, -5.97, 1.23])[0]\n",
      "\n",
      "#plotting\n",
      "fig0 = plt.figure(0)\n",
      "\n",
      "pts_ax = plt.plot(xarr,yarr,'b.')\n",
      "fig0.add_axes(plt.gca())\n",
      "\n",
      "xc = np.linspace(-9,-3,100)\n",
      "plt.plot(xc, Log10TransXVals(xc,*po),'b')\n",
      "fig0.add_axes(plt.gca())\n",
      "\n",
      "##############################\n",
      "#using NonLogXVals equation for fitting\n",
      "\n",
      "#fit the curve\n",
      "pof = curve_fit(NonLogXVals, xarr_10, yarr, [.49, 100.62, 10**-5.97, 1.23])[0]\n",
      "\n",
      "# plotting\n",
      "fig1 = plt.figure(1)\n",
      "## plot points\n",
      "plt.plot(xarr_10,yarr,'g.')\n",
      "fig1.add_axes(plt.gca())\n",
      "\n",
      "xf = np.logspace(-9,-3, 100)\n",
      "plt.plot(xf, NonLogXVals(xf,*pof),'g')\n",
      "plt.xscale('log')\n",
      "fig1.add_axes(plt.gca())\n",
      "\n",
      "\n",
      "################ \n",
      "# display to see that they produce the same results\n",
      "print \"Log10TransXVals() fitting results: \", po\n",
      "print \"NonLogXVals() fitting results: \", pof\n",
      "sys.stdout.flush()\n",
      "\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Log10TransXVals() fitting results:  [  -0.46229478  101.90361307   -5.92184104    1.10639417]\n",
        "NonLogXVals() fitting results:  [  1.01903602e+02  -4.62293730e-01   1.19717821e-06   1.10639466e+00]\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "fitting the data (Using Log10TransXVals function)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from matplotlib import rc\n",
      "import matplotlib as mpl\n",
      "fig_size = (10, 4)\n",
      "fig_dpi = 80*1.5\n",
      "plt.rc('font', family='sans-serif', size='12') \n",
      "mpl.rcParams['mathtext.default'] = 'regular'\n",
      "# to get sanserif, unitalicized math expression can also use\n",
      "# r'$\\mathrm{\\mathsf{'...'}}$' where latex math symbolic expression\n",
      "# replaces ellipses."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "'''\n",
      "mpl.text.Text(self, x=0, y=0, text=u'', color=None, verticalalignment=u'baseline', horizontalalignment=u'left', \n",
      "              multialignment=None, fontproperties=None, rotation=None, linespacing=None, rotation_mode=None, **kwargs)\n",
      "'''\n",
      "''"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "''"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# datapoints\n",
      "xlist = [-2.5,-2.,-1.5,-1.,-0.5, 0., .5, 1.] # exponents of x coords for both BzATP and ATP\n",
      "yb = [0.24886878,  0.29411765,  1.24434389,  3.25791855,  3.91402715] #BzATP y-vals\n",
      "ya = [0.29411765, 0.33936652,  0.54298643,  2.94117647,  4.77375566] #ATP y-vals\n",
      "\n",
      "#split up x coords (bc they only partially overlap for BzATP and ATP)\n",
      "xb = np.array(xlist[:5])\n",
      "xa = np.array(xlist[-5:])\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 125
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#function 1\n",
      "function = Log10TransXVals \n",
      "#fit the curves\n",
      "pb = curve_fit(function, xb, yb,[0.24886878,3.91402715,51*10**-3,1.23])[0]\n",
      "pa = curve_fit(function, xa, ya,[0.24886878,3.91402715,51*10**-3,1.23])[0]\n",
      "print pb, pa\n",
      "\n",
      "#make labels for fitted curves\n",
      "label_prep = lambda x: str(round(x,2))\n",
      "b_label = '$EC_{50}$ = %s $\\mu M$, h = %s'%(label_prep(10**pb[2]), label_prep(pb[3]))\n",
      "a_label = '$EC_{50}$ = %s $\\mu M$, h = %s'%(label_prep(10**pa[2]), label_prep(pa[3]))\n",
      "\n",
      "x_array = np.linspace(-3, 2, 70)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.20983406  3.98686534 -1.2966973   2.11169946] [ 0.31088598  4.90743233  0.45480966  2.79486493]\n"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#function 2\n",
      "function = NonLogXVals\n",
      "xb = 10**xb\n",
      "xa = 10**xa\n",
      "\n",
      "pb = curve_fit(function, xb, yb,[0.24886878,3.91402715,0.05050131653,1.23])[0]\n",
      "pa = curve_fit(function, xa, ya,[0.24886878,3.91402715,10**(2.88*10**-3),1.23])[0]\n",
      "print pb, pa\n",
      "\n",
      "#make labels for fitted curves\n",
      "label_prep = lambda x: str(round(x,2))\n",
      "b_label = '$EC_{50}$ = %s $\\mu M$, h = %s'%(label_prep(10**pb[2]), label_prep(pb[3]))\n",
      "a_label = '$EC_{50}$ = %s $\\mu M$, h = %s'%(label_prep(10**pa[2]), label_prep(pa[3]))\n",
      "\n",
      "x_array = np.logspace(-3, 2, 70)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 3.98686552  0.20983389  0.05050132  2.11169866] [ 0.31088598  4.90743233  2.849769   -2.79486496]\n"
       ]
      }
     ],
     "prompt_number": 126
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#plot the data points\n",
      "fig = plt.figure()\n",
      "if NonLogXVals == function:\n",
      "    plt.xlim(10**-3,10**2)\n",
      "    plt.xscale('log')\n",
      "else:\n",
      "    plt.xlim(-3,2)\n",
      "plt.plot(xb,yb,'bo',label='BzATP')\n",
      "plt.plot(xa,ya,'go', label='ATP')\n",
      "ax1 = fig.add_subplot(plt.gca())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 137
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "\n",
      "#plot fitted curves\n",
      "plt.plot(x_array, function(x_array,*pb), 'b',label=b_label)\n",
      "plt.plot(x_array, function(x_array,*pa), 'g',label=a_label)\n",
      "#plt.legend(loc=2,fontsize = 9.5)#'best')\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 138,
       "text": [
        "[<matplotlib.lines.Line2D at 0x156e1898>]"
       ]
      }
     ],
     "prompt_number": 138
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.ylim(0,5)\n",
      "# display and print\n",
      "print \"pb: \", pb, 10**pb[2],\n",
      "print \"pa: \", pa, 10**pa[2]\n",
      "sys.stdout.flush()\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "pb:  [ 3.98686552  0.20983389  0.05050132  2.11169866] 1.12331437195 pa:  [ 0.31088598  4.90743233  2.849769   -2.79486496] 707.569333424\n"
       ]
      }
     ],
     "prompt_number": 133
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ax2 = ax1.twiny()\n",
      "plt.xlim(10**-2.5,10)\n",
      "plt.xscale('log')\n",
      "\n",
      "plt.title('Fitting of Figure 1C in Yan et al. 2010*',y=1.059)#,fontdict={'fontsize':14})\n",
      "citation = 'Yan, Z., A. Khadra, S. Li, M. Tomic, A. Sherman, and S. S. Stojilkovic. (2010)\\\n",
      "Experimental Characterization and Mathematical Modeling of P2X7 Receptor Channel Gating. \\\n",
      "J Neurosci 30.42 : 14213-4224.'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 90
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Parasite axis"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from mpl_toolkits.axes_grid1 import host_subplot\n",
      "import mpl_toolkits.axisartist as AA\n",
      "import matplotlib.pyplot as plt\n",
      "\n",
      "host = host_subplot(111,axes_class=AA.Axes)\n",
      "plt.subplots_adjust(bottom=0.15)\n",
      "\n",
      "par1 = host.twiny()\n",
      "offset = -30\n",
      "new_fixed_axis = par1.get_grid_helper().new_fixed_axis\n",
      "par1.axis['bottom'] = new_fixed_axis(loc='bottom', axes=par1, offset=(0,offset))\n",
      "par1.axis['bottom'].toggle(all=True)\n",
      "\n",
      "host.set_xlim(-3,2)\n",
      "host.set_ylim(top=5,bottom=0)\n",
      "par1.set_xlim(10**-2.5,10)\n",
      "\n",
      "host.set_xlabel('Concentration (microMolars)')\n",
      "host.set_ylabel('Response')\n",
      "par1.set_xlabel('Concentration (log)')\n",
      "par1.set_xscale('log')\n",
      "\n",
      "p1, = host.plot(xb,yb,'bo',label='BzATP')\n",
      "p2, = host.plot(xa,ya,'go',label='ATP')\n",
      "p3, = host.plot(x_array, function(x_array,*pb), 'b',label=b_label)\n",
      "p4, = host.plot(x_array, function(x_array,*pa), 'g',label=a_label)\n",
      "host.set_xscale('log')\n",
      "\n",
      "host.legend(loc=2,fontsize = 9.5)#'best')\n",
      "\n",
      "plt.draw()\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "get x for one in terms of x for the other"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.lines as mlines\n",
      "import matplotlib.pyplot as plt\n",
      "\n",
      "# datapoints\n",
      "xlist = [-2.5,-2.,-1.5,-1.,-0.5, 0., .5, 1.] # exponents of x coords for both BzATP and ATP\n",
      "yb = [0.24886878,  0.29411765,  1.24434389,  3.25791855,  3.91402715] #BzATP y-vals\n",
      "ya = [0.29411765, 0.33936652,  0.54298643,  2.94117647,  4.77375566] #ATP y-vals\n",
      "\n",
      "#split up x coords (bc they only partially overlap for BzATP and ATP)\n",
      "xb = np.array(xlist[:5])\n",
      "xa = np.array(xlist[-5:])\n",
      "\n",
      "#function 2\n",
      "function = NonLogXVals\n",
      "xb = 10**xb\n",
      "xa = 10**xa\n",
      "\n",
      "pb = curve_fit(function, xb, yb,[0.24886878,3.91402715,0.05050131653,1.23])[0]\n",
      "pa = curve_fit(function, xa, ya,[0.24886878,3.91402715,10**(2.88*10**-3),1.23])[0]\n",
      "print pb, pa\n",
      "\n",
      "#make labels for fitted curves\n",
      "label_prep = lambda x: str(round(x,2))\n",
      "b_label = '$EC_{50}$ = %s , h = %s'%(label_prep(10**pb[2]), label_prep(pb[3]))\n",
      "a_label = '$EC_{50}$ = %s , h = %s'%(label_prep(10**pa[2]), label_prep(pa[3]))\n",
      "\n",
      "x_array = np.logspace(-3, 2, 70)\n",
      "\n",
      "\n",
      "#plot the data points\n",
      "fig = plt.figure()\n",
      "if NonLogXVals == function:\n",
      "    plt.xlim(10**-3,10**2)\n",
      "    plt.xscale('log')\n",
      "else:\n",
      "    plt.xlim(-3,2)\n",
      "plt.plot(xb,yb,'bo',label='BzATP')\n",
      "plt.plot(xa,ya,'go', label='ATP')\n",
      "ax1 = fig.add_subplot(plt.gca())\n",
      "\n",
      "\n",
      "#plot fitted curves\n",
      "plt.plot(x_array, function(x_array,*pb), 'b',label=b_label)\n",
      "plt.plot(x_array, function(x_array,*pa), 'g',label=a_label)\n",
      "#plt.legend(loc=2,fontsize = 9.5)#'best')\n",
      "\n",
      "\n",
      "atp_leg_key = mlines.Line2D([], [], color='g', marker='o', label='ATP (%s)'%a_label)\n",
      "bzatp_leg_key = mlines.Line2D([], [], color='b', marker='o', label='BzATP (%s)'%b_label)\n",
      "plt.legend(handles=[atp_leg_key,bzatp_leg_key], loc=2)\n",
      "\n",
      "\n",
      "#ATP\n",
      "'''\n",
      "a1=0.311\n",
      "b1=4.907\n",
      "ec1=2.850\n",
      "h1=-2.795\n",
      "'''\n",
      "a1,b1,ec1,h1 = pa\n",
      "\n",
      "#BzATP\n",
      "'''\n",
      "a2= 3.987\n",
      "b2= 0.210\n",
      "ec2=0.051\n",
      "h2=2.112\n",
      "'''\n",
      "a2,b2,ec2,h2 = pb\n",
      "\n",
      "def conversion(x2):\n",
      "    '''\n",
      "    num = a1 - b1\n",
      "    denom = b2 - b1 + ((a2-b2)/(1+(x2/ec2)**-h2))\n",
      "    x1 = ec1*((num/denom)-1)**-(1/h1)\n",
      "    '''\n",
      "    num = a1 - b1\n",
      "    denom = b2 - b1 + ((a2-b2)/(1+(x2/ec2)**-h2))\n",
      "    x1 = ec1*((num/denom)-1)**-(1/h1)\n",
      "    return x1\n",
      "\n",
      "\n",
      "sfreq = Slider(axfreq, '[BzATP]', 0.1, 30.0, valinit=2*10**-2)\n",
      "\n",
      "x_bzatp = 2*10**-2\n",
      "x_atp = conversion(x_bzatp)\n",
      "\n",
      "\n",
      "y_bzatp = function(x_bzatp,a2,b2,ec2,h2)\n",
      "y_atp = function(x_atp,a1,b1,ec1,h1)\n",
      "connect_line, = plt.plot([x_bzatp, x_atp],[y_bzatp, y_atp], 'k')\n",
      "b_point, = plt.plot(x_bzatp, y_bzatp,'b^')\n",
      "a_point, = plt.plot(x_atp, y_atp,'gv')\n",
      "\n",
      "\n",
      "\n",
      "plt.ylabel('I (nA)/40 (s)')\n",
      "plt.xlabel('Agonist Concentration (mM)')\n",
      "plt.title('Conversion from [BzATP] to [ATP] \\n(data from Yan et al. 2010*)')#,fontdict={'fontsize':14})\n",
      "\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 3.98686552  0.20983389  0.05050132  2.11169866] [ 0.31088598  4.90743233  2.849769   -2.79486496]\n"
       ]
      }
     ],
     "prompt_number": 169
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "----End of Curve Fitting Dose-Response Data for Logged and Unlogged X Values----"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, ax = plt.subplots()\n",
      "fig.canvas.draw_idle()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import pylab\n",
      "from scipy.optimize import curve_fit\n",
      " \n",
      "def sigmoid(x, x0, k):\n",
      "     y = 1 / (1 + np.exp(-k*(x-x0)))\n",
      "     return y\n",
      " \n",
      "xdata = np.array([0.0,   1.0,  3.0, 4.3, 7.0,   8.0,   8.5, 10.0, 12.0])\n",
      "ydata = np.array([0.01, 0.02, 0.04, 0.11, 0.43,  0.7, 0.89, 0.95, 0.99])\n",
      " \n",
      "popt, pcov = curve_fit(sigmoid, xdata, ydata)\n",
      "print popt\n",
      " \n",
      "x = np.linspace(-1, 15, 50)\n",
      "y = sigmoid(x, *popt)\n",
      " \n",
      "pylab.plot(xdata, ydata, 'o', label='data')\n",
      "pylab.plot(x,y, label='fit')\n",
      "pylab.ylim(0, 1.05)\n",
      "pylab.legend(loc='best')\n",
      "pylab.show()"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 7.14248671  1.12930053]\n"
       ]
      }
     ],
     "prompt_number": 44
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#fit gamma max\n",
      "from scipy.optimize import curve_fit\n",
      "from numpy import polyfit, set_printoptions\n",
      "set_printoptions(suppress=True)\n",
      "\n",
      "def line(x,a,b):\n",
      "    return a*x + b\n",
      "\n",
      "def poly2deg(x,c,b,a):\n",
      "    return c*(x**2) + b*(x) + a\n",
      "\n",
      "\n",
      "def poly3deg(x,d,c,b,a):\n",
      "    return d*(x**3) + c*(x**2) + b*(x) + a\n",
      "\n",
      "x = [-140,-120,-100,-80]\n",
      "y = [63,55,45,32]\n",
      "\n",
      "linAns = curve_fit(line, x, y) #[a b] = [-0.515 -7.9  ]\n",
      "print linAns\n",
      "print\n",
      "\n",
      "polyAns2 = polyfit(x,y,2)\n",
      "print polyAns2 #[ -0.003  -1.202 -44.15 ]\n",
      "\n",
      "polyAns3 = polyfit(x,y,3)\n",
      "print polyAns3 #[ -0.     -0.01   -1.942  -70.   ]\n",
      "\n",
      "v = linspace(-150,-60)\n",
      "plt.plot(x,y,'o-')\n",
      "plt.plot(v,line(v, -0.515, -7.9),'black')\n",
      "plt.plot(v,poly2deg(v,*polyAns2),'red')\n",
      "plt.plot(v,poly3deg(v,*polyAns3), 'green')\n",
      "\n",
      "plt.show()\n"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "\n",
      "\n",
      "zdelta = 0.34\n",
      "P0=0.11\n",
      "kboltzmann = 1.3806488e-23\n",
      "T =300\n",
      "nacl = 145\n",
      "gamma_max = \n",
      "\n",
      "ks =  lambda v:148*exp(-.21*v/(kboltzmann*T))\n",
      "gamma = lambda v: gamma_max/(1+(ks(v)/nacl))\n",
      "\n",
      "V = linspace(-10,10,.1)\n",
      "plt.plot(V,gamma(V))"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#trash\n",
      "def rearrangeStates(states):\n",
      "    '''\n",
      "    #states are organized first by time period, then by specific state type\n",
      "    #i.e.: if I run [ATP] = 3.2 for 20 seconds, then [ATP] = 0.0 for 10 more seconds\n",
      "    #states = [[[q1 vals],[q2 vals],...,[c4 vals]] for first 20 seconds,\n",
      "    #                    [[q1 vals],[q2 vals],...,[c4 vals]] for last 10 seconds]\n",
      "    #this function rearranges states so that they are organized only by state type \n",
      "    #i.e. qi from all time time periods get combined into one longer list of values\n",
      "    # returns [allQ1, allQ2, etc, allC1, etc]\n",
      "    '''\n",
      "    allStates = []\n",
      "    for states in zip(*states):\n",
      "        allStates.append(list(reduce(lambda state, element: itertools.chain(state,element),states)))\n",
      "    return allStates\n",
      "\n",
      "\n",
      "applyAg = 5\n",
      "washoutS = 45\n",
      "def agonist_concentration(t,atp):\n",
      "    if t >= applyAg and t < washoutS:\n",
      "        return atp\n",
      "    elif t < washoutS+40:\n",
      "        return atp+100\n",
      "    else:\n",
      "        return 0.0 "
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#construct function to return action for var\n",
      "def conditionals(conditions, actions):\n",
      "    if (len(conditions) != len(actions)) and not(len(conditions)+1 == len(actions)):\n",
      "        raise Exception(\n",
      "        (\"Conditions list has length {} and Actions lists has length {}.\"\n",
      "        \"Lists must either be the same length, or Actions list may contain\"\n",
      "        \"an action for when none of the conditions are met.\").format(len(conditions),len(actions))\n",
      "        )\n",
      "\n",
      "    def function(var):\n",
      "        ii=0\n",
      "        for cond,act in izip_longest(conditions,actions,fillvalue=lambda x:x==x):\n",
      "            ii += 1\n",
      "            #print ii,\" \",cond(var)\n",
      "            if cond(var):\n",
      "                return act(var)\n",
      "            else:\n",
      "                continue\n",
      "        #after iterating through all condition expressions provided by user, \n",
      "        #if none of the conditions are met and no 'safety' action has been specified\n",
      "        #at the end of the action list, raise an exception\n",
      "        raise Exception(\"Undefined case encountered. \" \n",
      "                        \"Provided conditions fail to cover all possible outcomes.\")\n",
      "    return function "
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#experimenting with ode\n",
      "\n",
      "tStart=0.\n",
      "tEnd=50.\n",
      "dose=32.*10**-4\n",
      "startConds=[0., 0., 0., 0., 1., 0., 0., 0.,Vs0, ns0, hs0]\n",
      "\n",
      "backend = 'dopri5'\n",
      "\n",
      "ode =  spi.ode(dfv).set_integrator(backend,nsteps=1)\n",
      "ode.set_initial_value(startConds,tStart).set_f_params(dose)\n",
      "# suppress Fortran-printed warning\n",
      "solver._integrator.iwork[2] = -1\n",
      "ts = []\n",
      "ys = []\n",
      "\n",
      "warnings.filterwarnings(\"ignore\", category=UserWarning)\n",
      "i=0\n",
      "while ode.t < tEnd:\n",
      "    try:\n",
      "        i += 1\n",
      "        if i%20000 == 0: print ode.t\n",
      "        ode.integrate(t1, step=True)\n",
      "        #ode.integrate(ode.t + dt)\n",
      "        #print \".\",\n",
      "    except Exception as e:\n",
      "        print e\n",
      "    ts.append(ode.t)\n",
      "    ys.append(ode.y)\n",
      "    if ode.successful():\n",
      "        print \"failure\"\n",
      "        q = raw_input()\n",
      "        if q == \"q\": break\n",
      "warnings.resetwarnings()\n",
      "print \"DONE\"\n",
      "\n",
      "ta = np.array(ts)\n",
      "ya = np.array(ys)\n",
      "print ya[0,]\n",
      "plt.figure(0)\n",
      "plt.plot(ta, ya.T[0],'r')\n",
      "plt.plot(ta, ya.T[1],'b')\n",
      "plt.plot(ta, ya.T[2],'orange')\n",
      "plt.plot(ta, ya.T[3],'yellow')\n",
      "plt.plot(ta, ya.T[4],'green')\n",
      "plt.plot(ta, ya.T[5],'purple')\n",
      "plt.plot(ta, ya.T[6],'black')\n",
      "plt.plot(ta, ya.T[7],'pink')\n",
      "#plt.show()\n",
      "plt.figure(1)\n",
      "plt.plot(ta, ya.T[-1],'g')\n",
      "plt.show()\n",
      "raw_input()"
     ],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {
      "run_control": {
       "read_only": false
      }
     },
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}