{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "####Glutamate Stimulation-will be used for normalization####\n",
      "import pickle\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "from neuron import h,gui\n",
      "\n",
      "####Sections and Connections####\n",
      "nofsections = 27 # 1soma + 25dendritic sections + 1terminating section \n",
      "nofspines =1     # number of spines is 1 throughout this simulation\n",
      "\n",
      "gc=[h.Section() for i in range(nofsections)]\n",
      "spineh=[h.Section() for i in range(nofspines)]\n",
      "spinen=[h.Section() for i in range(nofspines)]\n",
      "\n",
      "for i in range(nofsections-1):\n",
      "    gc[i+1].connect(gc[i],1,0)\n",
      "for i in range(nofspines):\n",
      "    spineh[i].connect(spinen[i],0,1)\n",
      "        \n",
      "j=0\n",
      "for i in range(nofspines): #spines start at 100um from soma, on 11th dendritic section\n",
      "    spinen[i].connect(gc[11+j],0.5,1)\n",
      "    j=j+1\n",
      "    \n",
      "####Morphology and other Parameters#### \n",
      "h.celsius=22 #temperature\n",
      "h.dt=0.025 #temporal resolution,ms\n",
      "\n",
      "for i in range(nofsections): \n",
      "    gc[i].L=10 #dendrite total length is 260 um\n",
      "    gc[i].nseg=3\n",
      "    gc[i].Ra=200\n",
      "    \n",
      "gc[0].diam=10 #soma size is 10umx10um, and there are two tapering regimes:\n",
      "for i in range(1,11):#1.tapering starts at 2.35um, ends at 1.7um , first 10 sections as in Ona-Jodar et al. Front Cell Neurosci 2017  \n",
      "    gc[i].diam=2.35-(i-1)*((2.35-1.7)/9)\n",
      "for i in range(11,27):#2.tapering starts at 1.7um ends at 1.2um, next 15 sections  \n",
      "    gc[i].diam=1.7-(i-10)*((1.7-1.2)/16)\n",
      "   \n",
      "for i in range(nofspines): \n",
      "    spineh[i].diam=1\n",
      "    spinen[i].diam=0.3 \n",
      "    spineh[i].L=1\n",
      "    spinen[i].L=2.5\n",
      "    spineh[i].nseg=3\n",
      "    spinen[i].nseg=3\n",
      "    spinen[i].Ra=4.9e3 #Ra is normalized as ohmcm.\n",
      "   \n",
      "presyn=h.Section() #with all other default specifications\n",
      "presyn.L=10\n",
      "presyn.diam=10\n",
      "\n",
      "####Settings for Ion channels and Synaptic Receptors and their Parameters#### \n",
      "for i in range(nofsections):\n",
      "    gc[i].insert('constant')#dummy current source\n",
      "    gc[i].insert('cadifusnpumpOGBenddif')#ca2+ and buffers diffusion,ca2+ pumps \n",
      "    gc[i].insert('nax')\n",
      "    gc[i].insert('kamt')\n",
      "    gc[i].insert('pas')\n",
      "    gc[i].g_pas=6e-4 \n",
      "    gc[i].e_pas=-85\n",
      "    gc[i].gbar_nax=0.5\n",
      "    gc[i].gbar_kamt=0.01\n",
      "    gc[i].cm=1\n",
      "    \n",
      "    gc[i].insert('canhem')#HVA Ca2+ channel \n",
      "    gc[i].insert('cathem')#T-type Ca2+ channel\n",
      "    gc[i].q10_cathem=3\n",
      "    gc[i].q10_canhem=3 \n",
      "    gc[i].a0m_cathem=0.055633 #to adjust the opening rate\n",
      "    gc[i].a0m_canhem=0.331432 #to adjust the opening rate   \n",
      "    gc[i].gcanbar_canhem=0.0005\n",
      "    gc[i].gcatbar_cathem=0.0003\n",
      "    \n",
      "    gc[i].TotalPump_cadifusnpumpOGBenddif=2e-11 \n",
      "    \n",
      "for i in range(nofspines):    \n",
      "    spineh[i].insert('constant')#dummy current source\n",
      "    spineh[i].insert('cadifusnpumpOGBenddif')#ca2+ and buffers diffusion,ca2+ pumps \n",
      "    spineh[i].insert('nax')\n",
      "    spineh[i].insert('kamt')\n",
      "    spineh[i].insert('pas')\n",
      "    spineh[i].gbar_nax=0.5 #for pharmacology in-silico,No Nav Channel case,set it to 0 \n",
      "    spineh[i].gbar_kamt=0.01\n",
      "    spineh[i].g_pas=2e-4\n",
      "    spineh[i].e_pas=-85\n",
      "    spineh[i].cm=1\n",
      "\n",
      "    spineh[i].insert('canhem')\n",
      "    spineh[i].insert('cathem')\n",
      "    spineh[i].q10_canhem=3\n",
      "    spineh[i].q10_cathem=3\n",
      "    spineh[i].a0m_cathem=0.055633 #to adjust the opening rate\n",
      "    spineh[i].a0m_canhem=0.331432 #to adjust the opening rate\n",
      "    spineh[i].gcatbar_cathem=0.00015 #for pharmacology in-silico,No T-type Ca2+ Channel case,set it to 0 \n",
      "    spineh[i].gcanbar_canhem=0.0004 #for pharmacology in-silico,No HVA Ca2+ Channel case,set it to 0 \n",
      "\n",
      "    spineh[i].TotalPump_cadifusnpumpOGBenddif=2.2e-11 \n",
      "    \n",
      "AMPARsyn=[h.AMPA5() for i in range(nofspines)]\n",
      "NMDARsyn=[h.NMDA5() for i in range(nofspines)]\n",
      "\n",
      "for i in range(nofspines): \n",
      "    AMPARsyn[i].loc(spineh[i](0.3))\n",
      "    AMPARsyn[i].gmax=2000\n",
      "    NMDARsyn[i].loc(spineh[i](0.7))\n",
      "    NMDARsyn[i].gmax=383 #for pharmacology in-silico,No NMDAR case,set it to 0.001 \n",
      "    NMDARsyn[i].gmax_ca=17 #for pharmacology in-silico,No NMDAR case,set it to 0.001 \n",
      "    ##NMDAR Setting##\n",
      "    NMDARsyn[i].Rb= 5e-3\n",
      "    NMDARsyn[i].Ru=12.9e-3\n",
      "    NMDARsyn[i].Rd=8.4e-3\n",
      "    NMDARsyn[i].Rr=6.8e-3\n",
      "    NMDARsyn[i].Ro=46.5e-3\n",
      "    NMDARsyn[i].Rc=73.8e-3 \n",
      "    ####\n",
      "    spinen[i].insert('cadifusnpumpOGBenddif')\n",
      "    spinen[i].TotalPump_cadifusnpumpOGBenddif=0 #there is no active mechanism on the neck\n",
      "\n",
      "####Setting Ca Dynamic Global Parameters####\n",
      "h.DCa_cadifusnpumpOGBenddif=0.6\n",
      "h.mg_NMDA5=1\n",
      "##endogenous buffer\n",
      "for i in range(nofsections):\n",
      "    gc[i].k1buf1_cadifusnpumpOGBenddif=1000\n",
      "    gc[i].k2buf1_cadifusnpumpOGBenddif=1\n",
      "    gc[i].TotalBuffer1_cadifusnpumpOGBenddif=0.12\n",
      "for i in range(nofspines):\n",
      "    spineh[i].k1buf1_cadifusnpumpOGBenddif=1000\n",
      "    spineh[i].k2buf1_cadifusnpumpOGBenddif=1\n",
      "    spineh[i].TotalBuffer1_cadifusnpumpOGBenddif=0.12\n",
      "##exogenous buffer    \n",
      "for i in range(nofsections):\n",
      "    gc[i].k1buf2_cadifusnpumpOGBenddif=1000\n",
      "    gc[i].k2buf2_cadifusnpumpOGBenddif=0.2\n",
      "    gc[i].TotalBuffer2_cadifusnpumpOGBenddif=0.1 #for \"no OGB case\", set this to 0\n",
      "for i in range(nofspines):\n",
      "    spineh[i].k1buf2_cadifusnpumpOGBenddif=1000\n",
      "    spineh[i].k2buf2_cadifusnpumpOGBenddif=0.2\n",
      "    spineh[i].TotalBuffer2_cadifusnpumpOGBenddif=0.1 #for \"no OGB case\", set this to 0\n",
      "\n",
      "####Mapping of Ca Concentration to fluorescence signal df/f \u2013 based on experimental data and simulations, see Figure 3C; not valid for \"no OGB case\"####\n",
      "def spine_fit(x):\n",
      "     y = 14390070 + (-49.1502 - 14390070)/(1 + (x/6632796000)**0.6715641)\n",
      "     return y\n",
      "    \n",
      "def dend_fit(x):  #not calculated and used based on experiment\n",
      "    y = 14390070 + (-49.1502 - 14390070)/(1 + (x/6632796000)**0.6715641)\n",
      "    return y\n",
      "\n",
      "####Simulation Readout####\n",
      "time_h = h.Vector()\n",
      "time_h.record(h._ref_t)\n",
      "vrec_gc=[h.Vector() for i in range(nofsections)] #gc[0] is the soma, gc[10] is the 1th parent dendrite\n",
      "vrec_spineh=[h.Vector() for i in range(nofspines)] \n",
      "icaspineh=[h.Vector() for i in range(nofspines)] #overall influx\n",
      "ccaspineh=[h.Vector() for i in range(nofspines)] #overall concentration\n",
      "icagc=[h.Vector() for i in range(nofsections)] #overall influx\n",
      "ccagc=[h.Vector() for i in range(nofsections)] #overall concentration\n",
      "\n",
      "icanmdr=[h.Vector() for i in range(nofspines)]\n",
      "icahvacc=[h.Vector() for i in range(nofspines)]\n",
      "icattype=[h.Vector() for i in range(nofspines)]\n",
      "\n",
      "##to check the buffering capacity- not used in the figures\n",
      "cbuf1spineh=[h.Vector() for i in range(nofspines)] \n",
      "ccabuf1spineh=[h.Vector() for i in range(nofspines)] \n",
      "cbuf2spineh=[h.Vector() for i in range(nofspines)] \n",
      "ccabuf2spineh=[h.Vector() for i in range(nofspines)]\n",
      "##\n",
      "for i in range(nofsections): \n",
      "    vrec_gc[i].record(gc[i](0.5)._ref_v)\n",
      "for i in range(nofspines): \n",
      "    vrec_spineh[i].record(spineh[i](0.5)._ref_v)    \n",
      "\n",
      "for i in range(nofspines): \n",
      "    icahvacc[i].record(spineh[i](0.5)._ref_ica_canhem)\n",
      "    icanmdr[i].record(NMDARsyn[i]._ref_ica)\n",
      "    icattype[i].record(spineh[i](0.5)._ref_ica_cathem)\n",
      "    \n",
      "    icaspineh[i].record(spineh[i](0.5)._ref_ica)\n",
      "    ccaspineh[i].record(spineh[i](0.5)._ref_caiav_cadifusnpumpOGBenddif)\n",
      "    \n",
      "    ##to check the buffering capacity- not used in the figures    \n",
      "    ccabuf1spineh[i].record(spineh[i](0.5)._ref_CaBufav1_cadifusnpumpOGBenddif)\n",
      "    ccabuf2spineh[i].record(spineh[i](0.5)._ref_CaBufav2_cadifusnpumpOGBenddif)\n",
      "    cbuf1spineh[i].record(spineh[i](0.5)._ref_Bufferav1_cadifusnpumpOGBenddif)\n",
      "    cbuf2spineh[i].record(spineh[i](0.5)._ref_Bufferav2_cadifusnpumpOGBenddif)\n",
      "    ##  \n",
      "    \n",
      "for i in range(nofsections):\n",
      "    icagc[i].record(gc[i](0.5)._ref_ica)\n",
      "    ccagc[i].record(gc[i](0.5)._ref_caiav_cadifusnpumpOGBenddif) \n",
      "    \n",
      "icanmdr_show=[np.array for i in range(nofspines)]\n",
      "icattype_show=[np.array for i in range(nofspines)]\n",
      "icahvacc_show=[np.array for i in range(nofspines)]\n",
      "ccaspineh_show=[np.array for i in range(nofspines)]\n",
      "icaspineh_show=[np.array for i in range(nofspines)]\n",
      "ccagc_show=[np.array for i in range(nofsections)]\n",
      "icagc_show=[np.array for i in range(nofsections)]    \n",
      "##to check the buffering capacity- not used in the figures \n",
      "cbuf1spineh_show=[np.array for i in range(nofsections)]\n",
      "cbuf2spineh_show=[np.array for i in range(nofsections)]\n",
      "ccabuf1spineh_show=[np.array for i in range(nofsections)]\n",
      "ccabuf2spineh_show=[np.array for i in range(nofsections)]\n",
      "##\n",
      "icanmdr_show=[np.array for i in range(nofspines)]\n",
      "icattype_show=[np.array for i in range(nofspines)]\n",
      "icahvacc_show=[np.array for i in range(nofspines)]\n",
      "y_dff_spineh=[np.array for i in range(nofspines)]\n",
      "y_dff_gc=[np.array for i in range(nofsections)]\n",
      "v_dend=[np.array for i in range(nofsections)]\n",
      "v_spineh=[np.array for i in range(nofspines)]\n",
      "\n",
      "vmspine_obj_glu=open(\"vmspine_0_record_glu\",\"w\")\n",
      "caspine_0_obj_glu=open(\"caspine_0_record_glu\",\"w\")\n",
      "dffspine_0_obj_glu=open(\"dffspine_0_record_glu\",\"w\")\n",
      "cagc_obj_glu=open(\"cagc_glu\",\"w\")\n",
      "dffgc_obj_glu=open(\"dffgc_glu\",\"w\")\n",
      "##to check the buffering capacity- not used in the figures \n",
      "ccabuf1spineh_obj_glu=open(\"ccabuf1spineh_glu\",\"w\")\n",
      "ccabuf2spineh_obj_glu=open(\"ccabuf2spineh_glu\",\"w\")\n",
      "cbuf1spineh_obj_glu=open(\"cbuf1spineh_glu\",\"w\")\n",
      "cbuf2spineh_obj_glu=open(\"cbuf2spineh_glu\",\"w\")                         \n",
      "##                          \n",
      "ittype_obj_glu=open(\"ittype_glu\",\"w\")\n",
      "ihvacc_obj_glu=open(\"ihvacc_glu\",\"w\")\n",
      "inmdr_obj_glu=open(\"inmdr_glu\",\"w\")\n",
      "\n",
      "t_c=[0]\n",
      "for tc in t_c:\n",
      "    \n",
      "####Setting Stimulation####\n",
      "#Glutamate\n",
      "    Rel=h.STEP_REL(0.75,presyn)\n",
      "    Rel.amplitude=1 \n",
      "    Rel.duration=3\n",
      "    Rel.release_time=120\n",
      "\n",
      "    for i in range(nofspines):\n",
      "        h.setpointer(Rel._ref_GLU,'C',AMPARsyn[i])\n",
      "        h.setpointer(Rel._ref_GLU,'C',NMDARsyn[i])\n",
      "    \n",
      "####Running the Simulation####\n",
      "    h.v_init=-85 #forced resting Vm for granule cells\n",
      "    h.init()\n",
      "\n",
      "    for l in range(nofsections):# dummy current source to compensate current caused by the forced Vm. \n",
      "        gc[l].ic_constant=-(gc[l].ina+gc[l].ik+gc[l].ica)\n",
      "    for l in range(nofspines):    \n",
      "        spineh[l].ic_constant=-(spineh[l].ina+spineh[l].ik+spineh[l].ica)\n",
      "\n",
      "    if h.cvode.active():\n",
      "        h.cvode.re_init()\n",
      "    else:\n",
      "        h.fcurrent()\n",
      "\n",
      "    h.tstop =350\n",
      "    h.run()\n",
      "    \n",
      "#### Vectors and conversion of units (um to nm)####\n",
      "    for i in range(nofspines):\n",
      "        v_spineh[i]=np.asarray(vrec_spineh[i])\n",
      "        ccaspineh_show[i]=1e6*np.asarray(ccaspineh[i])#converting to nM\n",
      "        icanmdr_show[i]=(1e2*np.asarray(icanmdr[i]))/3.14 #converting nA to mA/cm2, spine area is (piXe-8)cm2 \n",
      "        icahvacc_show[i]=np.asarray(icahvacc[i]) #mA/cm2\n",
      "        icattype_show[i]=np.asarray(icattype[i]) #mA/cm2\n",
      "        ##to check the buffering capacity- not used in the figures\n",
      "        ccabuf1spineh_show[i]=1e6*np.asarray(ccabuf1spineh[i])\n",
      "        ccabuf2spineh_show[i]=1e6*np.asarray(ccabuf2spineh[i])\n",
      "        cbuf1spineh_show[i]=1e6*np.asarray(cbuf1spineh[i])\n",
      "        cbuf2spineh_show[i]=1e6*np.asarray(cbuf2spineh[i])\n",
      "        ##\n",
      "    for i in range(nofsections):\n",
      "        v_dend[i]=np.asarray(vrec_gc[i])\n",
      "        ccagc_show[i]=1e6*np.asarray(ccagc[i])#converting to nM\n",
      "\n",
      "####Mapping of Ca Concentration to df/f####\n",
      "    for i in range(nofspines): \n",
      "        y_dff_spineh[i]=spine_fit(ccaspineh_show[i])\n",
      "    \n",
      "    for i in range(nofsections):    \n",
      "        y_dff_gc[i]=spine_fit(ccagc_show[i])    \n",
      "       \n",
      "    pickle.dump(v_spineh[0],vmspine_obj_glu)\n",
      "    pickle.dump(ccaspineh_show[0],caspine_0_obj_glu)\n",
      "    pickle.dump(y_dff_spineh[0],dffspine_0_obj_glu)\n",
      "    pickle.dump(ccagc_show[11],cagc_obj_glu)\n",
      "    pickle.dump(y_dff_gc[11],dffgc_obj_glu)\n",
      "    ##\n",
      "    pickle.dump(ccabuf1spineh_show[0],ccabuf1spineh_obj_glu)\n",
      "    pickle.dump(ccabuf2spineh_show[0],ccabuf2spineh_obj_glu)\n",
      "    pickle.dump(cbuf1spineh_show[0],cbuf1spineh_obj_glu)\n",
      "    pickle.dump(cbuf2spineh_show[0],cbuf2spineh_obj_glu)\n",
      "    ##\n",
      "    pickle.dump(icattype_show[0],ittype_obj_glu)\n",
      "    pickle.dump(icahvacc_show[0],ihvacc_obj_glu)\n",
      "    pickle.dump(icanmdr_show[0],inmdr_obj_glu)\n",
      "  \n",
      "caspine_0_obj_glu.close()\n",
      "dffspine_0_obj_glu.close()\n",
      "vmspine_obj_glu.close()\n",
      "cagc_obj_glu.close()\n",
      "dffgc_obj_glu.close()\n",
      "##\n",
      "ccabuf1spineh_obj_glu.close()\n",
      "ccabuf2spineh_obj_glu.close()\n",
      "cbuf1spineh_obj_glu.close()\n",
      "cbuf2spineh_obj_glu.close()\n",
      "##\n",
      "ittype_obj_glu.close()\n",
      "ihvacc_obj_glu.close()\n",
      "inmdr_obj_glu.close()\n",
      "\n",
      "print \"run the next block\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "end\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "####Current Clamp (global AP) Stimulation-will be used for normalization####\n",
      "\n",
      "import pickle\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "from neuron import h,gui\n",
      "\n",
      "####Sections and Connections####\n",
      "nofsections = 27 # 1soma + 25dendritic sections + 1terminating section \n",
      "nofspines =1     # number of spines is 1 throughout this simulation\n",
      "\n",
      "gc=[h.Section() for i in range(nofsections)]\n",
      "spineh=[h.Section() for i in range(nofspines)]\n",
      "spinen=[h.Section() for i in range(nofspines)]\n",
      "\n",
      "for i in range(nofsections-1):\n",
      "    gc[i+1].connect(gc[i],1,0)\n",
      "for i in range(nofspines):\n",
      "    spineh[i].connect(spinen[i],0,1)\n",
      "        \n",
      "j=0\n",
      "for i in range(nofspines): #spines start at 100um from soma, on 11th dendritic section\n",
      "    spinen[i].connect(gc[11+j],0.5,1)\n",
      "    j=j+1\n",
      "    \n",
      "####Morphology and other Parameters#### \n",
      "h.celsius=22 #temperature\n",
      "h.dt=0.025 #temporal resolution,ms\n",
      "\n",
      "for i in range(nofsections): \n",
      "    gc[i].L=10 #dendrite total length is 260 um\n",
      "    gc[i].nseg=3\n",
      "    gc[i].Ra=200\n",
      "    \n",
      "gc[0].diam=10 #soma size is 10umx10um, and there are two tapering regimes:\n",
      "for i in range(1,11):#1.tapering starts at 2.35um, ends at 1.7um , first 10 sections as in Ona-Jodar et al. Front Cell Neurosci 2017  \n",
      "    gc[i].diam=2.35-(i-1)*((2.35-1.7)/9)\n",
      "for i in range(11,27):#2.tapering starts at 1.7um ends at 1.2um, next 15 sections  \n",
      "    gc[i].diam=1.7-(i-10)*((1.7-1.2)/16)\n",
      "   \n",
      "for i in range(nofspines): \n",
      "    spineh[i].diam=1\n",
      "    spinen[i].diam=0.3 \n",
      "    spineh[i].L=1\n",
      "    spinen[i].L=2.5\n",
      "    spineh[i].nseg=3\n",
      "    spinen[i].nseg=3\n",
      "    spinen[i].Ra=4.9e3 #Ra is normalized as ohmcm.\n",
      "   \n",
      "presyn=h.Section() #with all other default specifications\n",
      "presyn.L=10\n",
      "presyn.diam=10\n",
      "\n",
      "####Settings for Ion channels and Synaptic Receptors and their Parameters#### \n",
      "for i in range(nofsections):\n",
      "    gc[i].insert('constant')#dummy current source\n",
      "    gc[i].insert('cadifusnpumpOGBenddif')#ca and buffers diffusion,ca pumps \n",
      "    gc[i].insert('nax')\n",
      "    gc[i].insert('kamt')\n",
      "    gc[i].insert('pas')\n",
      "    gc[i].g_pas=6e-4 \n",
      "    gc[i].e_pas=-85\n",
      "    gc[i].gbar_nax=0.5\n",
      "    gc[i].gbar_kamt=0.01\n",
      "    gc[i].cm=1\n",
      "    \n",
      "    gc[i].insert('canhem')#HVA Ca2+ channel \n",
      "    gc[i].insert('cathem')#T-type Ca2+ channel\n",
      "    gc[i].q10_cathem=3\n",
      "    gc[i].q10_canhem=3 \n",
      "    gc[i].a0m_cathem=0.055633 #to adjust the opening rate\n",
      "    gc[i].a0m_canhem=0.331432 #to adjust the opening rate   \n",
      "    gc[i].gcanbar_canhem=0.0005\n",
      "    gc[i].gcatbar_cathem=0.0003\n",
      "    \n",
      "    gc[i].TotalPump_cadifusnpumpOGBenddif=2e-11 \n",
      "    \n",
      "for i in range(nofspines):    \n",
      "    spineh[i].insert('constant')#dummy current source\n",
      "    spineh[i].insert('cadifusnpumpOGBenddif')#ca and buffers diffusion,ca pumps \n",
      "    spineh[i].insert('nax')\n",
      "    spineh[i].insert('kamt')\n",
      "    spineh[i].insert('pas')\n",
      "    spineh[i].gbar_nax=0.5 #for pharmacology in-silico,No Nav Channel case,set it to 0 \n",
      "    spineh[i].gbar_kamt=0.01\n",
      "    spineh[i].g_pas=2e-4\n",
      "    spineh[i].e_pas=-85\n",
      "    spineh[i].cm=1\n",
      "\n",
      "    spineh[i].insert('canhem')\n",
      "    spineh[i].insert('cathem')\n",
      "    spineh[i].q10_canhem=3\n",
      "    spineh[i].q10_cathem=3\n",
      "    spineh[i].a0m_cathem=0.055633 #to adjust the opening rate\n",
      "    spineh[i].a0m_canhem=0.331432 #to adjust the opening rate\n",
      "    spineh[i].gcatbar_cathem=0.00015#for pharmacology in-silico,No T-type Ca2+ Channel case,set it to 0 \n",
      "    spineh[i].gcanbar_canhem=0.0004 #for pharmacology in-silico,No HVA Ca2+ Channel case,set it to 0  \n",
      "\n",
      "    spineh[i].TotalPump_cadifusnpumpOGBenddif=2.2e-11 \n",
      "    \n",
      "AMPARsyn=[h.AMPA5() for i in range(nofspines)]\n",
      "NMDARsyn=[h.NMDA5() for i in range(nofspines)]\n",
      "\n",
      "for i in range(nofspines): \n",
      "    AMPARsyn[i].loc(spineh[i](0.3))\n",
      "    AMPARsyn[i].gmax=2000\n",
      "    NMDARsyn[i].loc(spineh[i](0.7))\n",
      "    NMDARsyn[i].gmax=383 #for pharmacology in-silico,No NMDAR case,set it to 0.001\n",
      "    NMDARsyn[i].gmax_ca=17 #for pharmacology in-silico,No NMDAR case,set it to 0.001 \n",
      "    ##NMDAR Setting##\n",
      "    NMDARsyn[i].Rb= 5e-3\n",
      "    NMDARsyn[i].Ru=12.9e-3\n",
      "    NMDARsyn[i].Rd=8.4e-3\n",
      "    NMDARsyn[i].Rr=6.8e-3\n",
      "    NMDARsyn[i].Ro=46.5e-3\n",
      "    NMDARsyn[i].Rc=73.8e-3 \n",
      "    ####\n",
      "    spinen[i].insert('cadifusnpumpOGBenddif')\n",
      "    spinen[i].TotalPump_cadifusnpumpOGBenddif=0 #there is no active mechanism on the neck\n",
      "\n",
      "####Setting Ca Dynamic Global Parameters####\n",
      "h.DCa_cadifusnpumpOGBenddif=0.6\n",
      "h.mg_NMDA5=1\n",
      "##endogenous buffer\n",
      "for i in range(nofsections):\n",
      "    gc[i].k1buf1_cadifusnpumpOGBenddif=1000\n",
      "    gc[i].k2buf1_cadifusnpumpOGBenddif=1\n",
      "    gc[i].TotalBuffer1_cadifusnpumpOGBenddif=0.12\n",
      "for i in range(nofspines):\n",
      "    spineh[i].k1buf1_cadifusnpumpOGBenddif=1000\n",
      "    spineh[i].k2buf1_cadifusnpumpOGBenddif=1\n",
      "    spineh[i].TotalBuffer1_cadifusnpumpOGBenddif=0.12\n",
      "##exogenous buffer    \n",
      "for i in range(nofsections):\n",
      "    gc[i].k1buf2_cadifusnpumpOGBenddif=1000\n",
      "    gc[i].k2buf2_cadifusnpumpOGBenddif=0.2\n",
      "    gc[i].TotalBuffer2_cadifusnpumpOGBenddif=0.1 #for \"no OGB case\", set this to 0\n",
      "for i in range(nofspines):\n",
      "    spineh[i].k1buf2_cadifusnpumpOGBenddif=1000\n",
      "    spineh[i].k2buf2_cadifusnpumpOGBenddif=0.2\n",
      "    spineh[i].TotalBuffer2_cadifusnpumpOGBenddif=0.1 #for \"no OGB case\", set this to 0\n",
      "\n",
      "####Mapping of Ca Concentration to fluorescence signal df/f \u2013 based on experimental data and simulations, see Figure 3C; not valid for \"no OGB case\"####\n",
      "def spine_fit(x):\n",
      "     y = 14390070 + (-49.1502 - 14390070)/(1 + (x/6632796000)**0.6715641)\n",
      "     return y\n",
      "    \n",
      "def dend_fit(x):  #not calculated and used based on experiment\n",
      "    y = 14390070 + (-49.1502 - 14390070)/(1 + (x/6632796000)**0.6715641)\n",
      "    return y\n",
      "\n",
      "####Simulation Readout####\n",
      "time_h = h.Vector()\n",
      "time_h.record(h._ref_t)\n",
      "vrec_gc=[h.Vector() for i in range(nofsections)] #gc[0] is the soma, gc[10] is the 1th parent dendrite\n",
      "vrec_spineh=[h.Vector() for i in range(nofspines)] \n",
      "icaspineh=[h.Vector() for i in range(nofspines)] #overall influx\n",
      "ccaspineh=[h.Vector() for i in range(nofspines)] #overall concentration\n",
      "icagc=[h.Vector() for i in range(nofsections)] #overall influx\n",
      "ccagc=[h.Vector() for i in range(nofsections)] #overall concentration\n",
      "\n",
      "icahvacc=[h.Vector() for i in range(nofspines)]\n",
      "icattype=[h.Vector() for i in range(nofspines)]\n",
      "\n",
      "##to check the buffering capacity- not used in the figures\n",
      "cbuf1spineh=[h.Vector() for i in range(nofspines)] \n",
      "ccabuf1spineh=[h.Vector() for i in range(nofspines)] \n",
      "cbuf2spineh=[h.Vector() for i in range(nofspines)] \n",
      "ccabuf2spineh=[h.Vector() for i in range(nofspines)]\n",
      "##\n",
      "for i in range(nofsections): \n",
      "    vrec_gc[i].record(gc[i](0.5)._ref_v)\n",
      "for i in range(nofspines): \n",
      "    vrec_spineh[i].record(spineh[i](0.5)._ref_v)    \n",
      "\n",
      "for i in range(nofspines): \n",
      "    icahvacc[i].record(spineh[i](0.5)._ref_ica_canhem)\n",
      "    icattype[i].record(spineh[i](0.5)._ref_ica_cathem)\n",
      "    \n",
      "    icaspineh[i].record(spineh[i](0.5)._ref_ica)\n",
      "    ccaspineh[i].record(spineh[i](0.5)._ref_caiav_cadifusnpumpOGBenddif)\n",
      "    \n",
      "    ##to check the buffering capacity- not used in the figures    \n",
      "    ccabuf1spineh[i].record(spineh[i](0.5)._ref_CaBufav1_cadifusnpumpOGBenddif)\n",
      "    ccabuf2spineh[i].record(spineh[i](0.5)._ref_CaBufav2_cadifusnpumpOGBenddif)\n",
      "    cbuf1spineh[i].record(spineh[i](0.5)._ref_Bufferav1_cadifusnpumpOGBenddif)\n",
      "    cbuf2spineh[i].record(spineh[i](0.5)._ref_Bufferav2_cadifusnpumpOGBenddif)\n",
      "    ##  \n",
      "    \n",
      "for i in range(nofsections):\n",
      "    icagc[i].record(gc[i](0.5)._ref_ica)\n",
      "    ccagc[i].record(gc[i](0.5)._ref_caiav_cadifusnpumpOGBenddif) \n",
      "    \n",
      "icattype_show=[np.array for i in range(nofspines)]\n",
      "icahvacc_show=[np.array for i in range(nofspines)]\n",
      "ccaspineh_show=[np.array for i in range(nofspines)]\n",
      "icaspineh_show=[np.array for i in range(nofspines)]\n",
      "ccagc_show=[np.array for i in range(nofsections)]\n",
      "icagc_show=[np.array for i in range(nofsections)]    \n",
      "##to check the buffering capacity- not used in the figures \n",
      "cbuf1spineh_show=[np.array for i in range(nofsections)]\n",
      "cbuf2spineh_show=[np.array for i in range(nofsections)]\n",
      "ccabuf1spineh_show=[np.array for i in range(nofsections)]\n",
      "ccabuf2spineh_show=[np.array for i in range(nofsections)]\n",
      "##\n",
      "icattype_show=[np.array for i in range(nofspines)]\n",
      "icahvacc_show=[np.array for i in range(nofspines)]\n",
      "y_dff_spineh=[np.array for i in range(nofspines)]\n",
      "y_dff_gc=[np.array for i in range(nofsections)]\n",
      "v_dend=[np.array for i in range(nofsections)]\n",
      "v_spineh=[np.array for i in range(nofspines)]\n",
      "\n",
      "vmspine_obj_ap=open(\"vmspine_0_record_ap\",\"w\")\n",
      "caspine_0_obj_ap=open(\"caspine_0_record_ap\",\"w\")\n",
      "dffspine_0_obj_ap=open(\"dffspine_0_record_ap\",\"w\")\n",
      "cagc_obj_ap=open(\"cagc_ap\",\"w\")\n",
      "dffgc_obj_ap=open(\"dffgc_ap\",\"w\")\n",
      "##to check the buffering capacity- not used in the figures \n",
      "ccabuf1spineh_obj_ap=open(\"ccabuf1spineh_ap\",\"w\")\n",
      "ccabuf2spineh_obj_ap=open(\"ccabuf2spineh_ap\",\"w\")\n",
      "cbuf1spineh_obj_ap=open(\"cbuf1spineh_ap\",\"w\")\n",
      "cbuf2spineh_obj_ap=open(\"cbuf2spineh_ap\",\"w\")                         \n",
      "##                          \n",
      "ittype_obj_ap=open(\"ittype_ap\",\"w\")\n",
      "ihvacc_obj_ap=open(\"ihvacc_ap\",\"w\")\n",
      "\n",
      "t_c=[-100,-80,-60,-40,-20,-12,-10,-8,-7,-6,-5,-4,-3,-2,0,2,3,4,5,6,7,8,10,12,20,40,60,80,100]\n",
      "for tc in t_c:\n",
      "####Setting Stimulation####\n",
      "#Current Clamp (Global AP)\n",
      "    APstim1=h.IClamp(0.5,sec=gc[0])\n",
      "    APstim1.delay=120+tc #mS\n",
      "    APstim1.dur=3 #mS\n",
      "    APstim1.amp=1 #nA\n",
      "#Glutamate-not applicable in this block, just mentioned to avoid pointer problem.amplitude=0 \n",
      "    Rel=h.STEP_REL(0.75,presyn)\n",
      "    Rel.amplitude=0\n",
      "    Rel.duration=3\n",
      "    Rel.release_time=120\n",
      "\n",
      "    for i in range(nofspines):\n",
      "        h.setpointer(Rel._ref_GLU,'C',AMPARsyn[i])\n",
      "        h.setpointer(Rel._ref_GLU,'C',NMDARsyn[i])\n",
      "    \n",
      "####Running the Simulation####\n",
      "    h.v_init=-85 #forced resting Vm for granule cells\n",
      "    h.init()\n",
      "\n",
      "    for l in range(nofsections):# dummy current source to compensate current caused by the forced Vm. \n",
      "        gc[l].ic_constant=-(gc[l].ina+gc[l].ik+gc[l].ica)\n",
      "    for l in range(nofspines):    \n",
      "        spineh[l].ic_constant=-(spineh[l].ina+spineh[l].ik+spineh[l].ica)\n",
      "\n",
      "    if h.cvode.active():\n",
      "        h.cvode.re_init()\n",
      "    else:\n",
      "        h.fcurrent()\n",
      "\n",
      "    h.tstop =350\n",
      "    h.run()\n",
      "    \n",
      "#### Vectors and conversion of units (um to nm)####\n",
      "    for i in range(nofspines):\n",
      "        v_spineh[i]=np.asarray(vrec_spineh[i])\n",
      "        ccaspineh_show[i]=1e6*np.asarray(ccaspineh[i])#converting to nM\n",
      "        icahvacc_show[i]=np.asarray(icahvacc[i]) #mA/cm2\n",
      "        icattype_show[i]=np.asarray(icattype[i]) #mA/cm2\n",
      "        ##to check the buffering capacity- not used in the figures\n",
      "        ccabuf1spineh_show[i]=1e6*np.asarray(ccabuf1spineh[i])\n",
      "        ccabuf2spineh_show[i]=1e6*np.asarray(ccabuf2spineh[i])\n",
      "        cbuf1spineh_show[i]=1e6*np.asarray(cbuf1spineh[i])\n",
      "        cbuf2spineh_show[i]=1e6*np.asarray(cbuf2spineh[i])\n",
      "        ##\n",
      "    for i in range(nofsections):\n",
      "        v_dend[i]=np.asarray(vrec_gc[i])\n",
      "        ccagc_show[i]=1e6*np.asarray(ccagc[i])#converting to nM\n",
      "\n",
      "####Mapping of Ca Concentration to df/f####\n",
      "    for i in range(nofspines): \n",
      "        y_dff_spineh[i]=spine_fit(ccaspineh_show[i])\n",
      "    \n",
      "    for i in range(nofsections):    \n",
      "        y_dff_gc[i]=spine_fit(ccagc_show[i])    \n",
      "       \n",
      "    pickle.dump(v_spineh[0],vmspine_obj_ap)\n",
      "    pickle.dump(ccaspineh_show[0],caspine_0_obj_ap)\n",
      "    pickle.dump(y_dff_spineh[0],dffspine_0_obj_ap)\n",
      "    pickle.dump(ccagc_show[11],cagc_obj_ap)\n",
      "    pickle.dump(y_dff_gc[11],dffgc_obj_ap)\n",
      "    ##\n",
      "    pickle.dump(ccabuf1spineh_show[0],ccabuf1spineh_obj_ap)\n",
      "    pickle.dump(ccabuf2spineh_show[0],ccabuf2spineh_obj_ap)\n",
      "    pickle.dump(cbuf1spineh_show[0],cbuf1spineh_obj_ap)\n",
      "    pickle.dump(cbuf2spineh_show[0],cbuf2spineh_obj_ap)\n",
      "    ##\n",
      "    pickle.dump(icattype_show[0],ittype_obj_ap)\n",
      "    pickle.dump(icahvacc_show[0],ihvacc_obj_ap)\n",
      "  \n",
      "caspine_0_obj_ap.close()\n",
      "dffspine_0_obj_ap.close()\n",
      "vmspine_obj_ap.close()\n",
      "cagc_obj_ap.close()\n",
      "dffgc_obj_ap.close()\n",
      "##\n",
      "ccabuf1spineh_obj_ap.close()\n",
      "ccabuf2spineh_obj_ap.close()\n",
      "cbuf1spineh_obj_ap.close()\n",
      "cbuf2spineh_obj_ap.close()\n",
      "##\n",
      "ittype_obj_ap.close()\n",
      "ihvacc_obj_ap.close()\n",
      "\n",
      "print \"run the next block\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "run the next block\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "####Coincidance Case-Glutamate and Global AP####\n",
      "\n",
      "import pickle\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "from neuron import h,gui\n",
      "\n",
      "####Sections and Connections####\n",
      "nofsections = 27 # 1soma + 25dendritic sections + 1terminating section \n",
      "nofspines =1     # number of spines is 1 throughout this simulation\n",
      "\n",
      "gc=[h.Section() for i in range(nofsections)]\n",
      "spineh=[h.Section() for i in range(nofspines)]\n",
      "spinen=[h.Section() for i in range(nofspines)]\n",
      "\n",
      "for i in range(nofsections-1):\n",
      "    gc[i+1].connect(gc[i],1,0)\n",
      "for i in range(nofspines):\n",
      "    spineh[i].connect(spinen[i],0,1)\n",
      "        \n",
      "j=0\n",
      "for i in range(nofspines): #spines start at 100um from soma, on 11th dendritic section\n",
      "    spinen[i].connect(gc[11+j],0.5,1)\n",
      "    j=j+1\n",
      "    \n",
      "####Morphology and other Parameters#### \n",
      "h.celsius=22 #temperature\n",
      "h.dt=0.025 #temporal resolution,ms\n",
      "\n",
      "for i in range(nofsections): \n",
      "    gc[i].L=10 #dendrite total length is 260 um\n",
      "    gc[i].nseg=3\n",
      "    gc[i].Ra=200\n",
      "    \n",
      "gc[0].diam=10 #soma size is 10umx10um, and there are two tapering regimes:\n",
      "for i in range(1,11):#1.tapering starts at 2.35um, ends at 1.7um , first 10 sections as in Ona-Jodar et al. Front Cell Neurosci 2017  \n",
      "    gc[i].diam=2.35-(i-1)*((2.35-1.7)/9)\n",
      "for i in range(11,27):#2.tapering starts at 1.7um ends at 1.2um, next 15 sections  \n",
      "    gc[i].diam=1.7-(i-10)*((1.7-1.2)/16)\n",
      "   \n",
      "for i in range(nofspines): \n",
      "    spineh[i].diam=1\n",
      "    spinen[i].diam=0.3 \n",
      "    spineh[i].L=1\n",
      "    spinen[i].L=2.5\n",
      "    spineh[i].nseg=3\n",
      "    spinen[i].nseg=3\n",
      "    spinen[i].Ra=4.9e3 #Ra is normalized as ohmcm.\n",
      "   \n",
      "presyn=h.Section() #with all other default specifications\n",
      "presyn.L=10\n",
      "presyn.diam=10\n",
      "\n",
      "####Settings for Ion channels and Synaptic Receptors and their Parameters#### \n",
      "for i in range(nofsections):\n",
      "    gc[i].insert('constant')#dummy current source\n",
      "    gc[i].insert('cadifusnpumpOGBenddif')#ca and buffers diffusion,ca pumps \n",
      "    gc[i].insert('nax')\n",
      "    gc[i].insert('kamt')\n",
      "    gc[i].insert('pas')\n",
      "    gc[i].g_pas=6e-4 \n",
      "    gc[i].e_pas=-85\n",
      "    gc[i].gbar_nax=0.5\n",
      "    gc[i].gbar_kamt=0.01\n",
      "    gc[i].cm=1\n",
      "    \n",
      "    gc[i].insert('canhem')#HVA Ca2+ channel \n",
      "    gc[i].insert('cathem')#T-type Ca2+ channel\n",
      "    gc[i].q10_cathem=3\n",
      "    gc[i].q10_canhem=3 \n",
      "    gc[i].a0m_cathem=0.055633 #to adjust the opening rate\n",
      "    gc[i].a0m_canhem=0.331432 #to adjust the opening rate   \n",
      "    gc[i].gcanbar_canhem=0.0005\n",
      "    gc[i].gcatbar_cathem=0.0003\n",
      "    \n",
      "    gc[i].TotalPump_cadifusnpumpOGBenddif=2e-11 \n",
      "    \n",
      "for i in range(nofspines):    \n",
      "    spineh[i].insert('constant')#dummy current source\n",
      "    spineh[i].insert('cadifusnpumpOGBenddif')#ca and buffers diffusion,ca pumps \n",
      "    spineh[i].insert('nax')\n",
      "    spineh[i].insert('kamt')\n",
      "    spineh[i].insert('pas')\n",
      "    spineh[i].gbar_nax=0.5 #for pharmacology in-silico,No Nav Channel case,set it to 0 \n",
      "    spineh[i].gbar_kamt=0.01\n",
      "    spineh[i].g_pas=2e-4\n",
      "    spineh[i].e_pas=-85\n",
      "    spineh[i].cm=1\n",
      "\n",
      "    spineh[i].insert('canhem')\n",
      "    spineh[i].insert('cathem')\n",
      "    spineh[i].q10_canhem=3\n",
      "    spineh[i].q10_cathem=3\n",
      "    spineh[i].a0m_cathem=0.055633 #to adjust the opening rate\n",
      "    spineh[i].a0m_canhem=0.331432 #to adjust the opening rate\n",
      "    spineh[i].gcatbar_cathem=0.00015 #for pharmacology in-silico,No T-type Ca2+ Channel case,set it to 0\n",
      "    spineh[i].gcanbar_canhem=0.0004 #for pharmacology in-silico,No HVA Ca2+ Channel case,set it to 0 \n",
      "\n",
      "    spineh[i].TotalPump_cadifusnpumpOGBenddif=2.2e-11 \n",
      "    \n",
      "AMPARsyn=[h.AMPA5() for i in range(nofspines)]\n",
      "NMDARsyn=[h.NMDA5() for i in range(nofspines)]\n",
      "\n",
      "for i in range(nofspines): \n",
      "    AMPARsyn[i].loc(spineh[i](0.3))\n",
      "    AMPARsyn[i].gmax=2000\n",
      "    NMDARsyn[i].loc(spineh[i](0.7))\n",
      "    NMDARsyn[i].gmax=383 #for pharmacology in-silico,No NMDAR case,set it to 0.001\n",
      "    NMDARsyn[i].gmax_ca=17#for pharmacology in-silico,No NMDAR case,set it to 0.001 \n",
      "    ##NMDAR Setting##\n",
      "    NMDARsyn[i].Rb= 5e-3\n",
      "    NMDARsyn[i].Ru=12.9e-3\n",
      "    NMDARsyn[i].Rd=8.4e-3\n",
      "    NMDARsyn[i].Rr=6.8e-3\n",
      "    NMDARsyn[i].Ro=46.5e-3\n",
      "    NMDARsyn[i].Rc=73.8e-3 \n",
      "    ####\n",
      "    spinen[i].insert('cadifusnpumpOGBenddif')\n",
      "    spinen[i].TotalPump_cadifusnpumpOGBenddif=0 #there is no active mechanism on the neck\n",
      "\n",
      "####Setting Ca Dynamic Global Parameters####\n",
      "h.DCa_cadifusnpumpOGBenddif=0.6\n",
      "h.mg_NMDA5=1\n",
      "##endogenous buffer\n",
      "for i in range(nofsections):\n",
      "    gc[i].k1buf1_cadifusnpumpOGBenddif=1000\n",
      "    gc[i].k2buf1_cadifusnpumpOGBenddif=1\n",
      "    gc[i].TotalBuffer1_cadifusnpumpOGBenddif=0.12\n",
      "for i in range(nofspines):\n",
      "    spineh[i].k1buf1_cadifusnpumpOGBenddif=1000\n",
      "    spineh[i].k2buf1_cadifusnpumpOGBenddif=1\n",
      "    spineh[i].TotalBuffer1_cadifusnpumpOGBenddif=0.12\n",
      "##exogenous buffer    \n",
      "for i in range(nofsections):\n",
      "    gc[i].k1buf2_cadifusnpumpOGBenddif=1000\n",
      "    gc[i].k2buf2_cadifusnpumpOGBenddif=0.2\n",
      "    gc[i].TotalBuffer2_cadifusnpumpOGBenddif=0.1 #for \"no OGB case\", set this to 0\n",
      "for i in range(nofspines):\n",
      "    spineh[i].k1buf2_cadifusnpumpOGBenddif=1000\n",
      "    spineh[i].k2buf2_cadifusnpumpOGBenddif=0.2\n",
      "    spineh[i].TotalBuffer2_cadifusnpumpOGBenddif=0.1 #for \"no OGB case\", set this to 0\n",
      "\n",
      "####Mapping of Ca Concentration to fluorescence signal df/f \u2013 based on experimental data and simulations, see Figure 3C; not valid for \"no OGB case\"####\n",
      "def spine_fit(x):\n",
      "     y = 14390070 + (-49.1502 - 14390070)/(1 + (x/6632796000)**0.6715641)\n",
      "     return y\n",
      "    \n",
      "def dend_fit(x):  #not calculated and used based on experiment\n",
      "    y = 14390070 + (-49.1502 - 14390070)/(1 + (x/6632796000)**0.6715641)\n",
      "    return y\n",
      "\n",
      "####Simulation Readout####\n",
      "time_h = h.Vector()\n",
      "time_h.record(h._ref_t)\n",
      "vrec_gc=[h.Vector() for i in range(nofsections)] #gc[0] is the soma, gc[10] is the 1th parent dendrite\n",
      "vrec_spineh=[h.Vector() for i in range(nofspines)] \n",
      "icaspineh=[h.Vector() for i in range(nofspines)] #overall influx\n",
      "ccaspineh=[h.Vector() for i in range(nofspines)] #overall concentration\n",
      "icagc=[h.Vector() for i in range(nofsections)] #overall influx\n",
      "ccagc=[h.Vector() for i in range(nofsections)] #overall concentration\n",
      "\n",
      "icahvacc=[h.Vector() for i in range(nofspines)]\n",
      "icattype=[h.Vector() for i in range(nofspines)]\n",
      "\n",
      "##to check the buffering capacity- not used in the figures\n",
      "cbuf1spineh=[h.Vector() for i in range(nofspines)] \n",
      "ccabuf1spineh=[h.Vector() for i in range(nofspines)] \n",
      "cbuf2spineh=[h.Vector() for i in range(nofspines)] \n",
      "ccabuf2spineh=[h.Vector() for i in range(nofspines)]\n",
      "##\n",
      "for i in range(nofsections): \n",
      "    vrec_gc[i].record(gc[i](0.5)._ref_v)\n",
      "for i in range(nofspines): \n",
      "    vrec_spineh[i].record(spineh[i](0.5)._ref_v)    \n",
      "\n",
      "for i in range(nofspines): \n",
      "    icahvacc[i].record(spineh[i](0.5)._ref_ica_canhem)\n",
      "    icattype[i].record(spineh[i](0.5)._ref_ica_cathem)\n",
      "    \n",
      "    icaspineh[i].record(spineh[i](0.5)._ref_ica)\n",
      "    ccaspineh[i].record(spineh[i](0.5)._ref_caiav_cadifusnpumpOGBenddif)\n",
      "    \n",
      "    ##to check the buffering capacity- not used in the figures    \n",
      "    ccabuf1spineh[i].record(spineh[i](0.5)._ref_CaBufav1_cadifusnpumpOGBenddif)\n",
      "    ccabuf2spineh[i].record(spineh[i](0.5)._ref_CaBufav2_cadifusnpumpOGBenddif)\n",
      "    cbuf1spineh[i].record(spineh[i](0.5)._ref_Bufferav1_cadifusnpumpOGBenddif)\n",
      "    cbuf2spineh[i].record(spineh[i](0.5)._ref_Bufferav2_cadifusnpumpOGBenddif)\n",
      "    ##  \n",
      "    \n",
      "for i in range(nofsections):\n",
      "    icagc[i].record(gc[i](0.5)._ref_ica)\n",
      "    ccagc[i].record(gc[i](0.5)._ref_caiav_cadifusnpumpOGBenddif) \n",
      "    \n",
      "icattype_show=[np.array for i in range(nofspines)]\n",
      "icahvacc_show=[np.array for i in range(nofspines)]\n",
      "ccaspineh_show=[np.array for i in range(nofspines)]\n",
      "icaspineh_show=[np.array for i in range(nofspines)]\n",
      "ccagc_show=[np.array for i in range(nofsections)]\n",
      "icagc_show=[np.array for i in range(nofsections)]    \n",
      "##to check the buffering capacity- not used in the figures \n",
      "cbuf1spineh_show=[np.array for i in range(nofsections)]\n",
      "cbuf2spineh_show=[np.array for i in range(nofsections)]\n",
      "ccabuf1spineh_show=[np.array for i in range(nofsections)]\n",
      "ccabuf2spineh_show=[np.array for i in range(nofsections)]\n",
      "##\n",
      "icattype_show=[np.array for i in range(nofspines)]\n",
      "icahvacc_show=[np.array for i in range(nofspines)]\n",
      "y_dff_spineh=[np.array for i in range(nofspines)]\n",
      "y_dff_gc=[np.array for i in range(nofsections)]\n",
      "v_dend=[np.array for i in range(nofsections)]\n",
      "v_spineh=[np.array for i in range(nofspines)]\n",
      "\n",
      "vmspine_obj=open(\"vmspine_0_record\",\"w\")\n",
      "caspine_0_obj=open(\"caspine_0_record\",\"w\")\n",
      "dffspine_0_obj=open(\"dffspine_0_record\",\"w\")\n",
      "cagc_obj=open(\"cagc\",\"w\")\n",
      "dffgc_obj=open(\"dffgc\",\"w\")\n",
      "##to check the buffering capacity- not used in the figures \n",
      "ccabuf1spineh_obj=open(\"ccabuf1spineh\",\"w\")\n",
      "ccabuf2spineh_obj=open(\"ccabuf2spineh\",\"w\")\n",
      "cbuf1spineh_obj=open(\"cbuf1spineh\",\"w\")\n",
      "cbuf2spineh_obj=open(\"cbuf2spineh\",\"w\")                         \n",
      "##                          \n",
      "ittype_obj=open(\"ittype\",\"w\")\n",
      "ihvacc_obj=open(\"ihvacc\",\"w\")\n",
      "\n",
      "t_c=[-100,-80,-60,-40,-20,-12,-10,-8,-7,-6,-5,-4,-3,-2,0,2,3,4,5,6,7,8,10,12,20,40,60,80,100]\n",
      "for tc in t_c:\n",
      "####Setting Stimulation####\n",
      "#Current Clamp (Global AP)\n",
      "    APstim1=h.IClamp(0.5,sec=gc[0])\n",
      "    APstim1.delay=120+tc #mS\n",
      "    APstim1.dur=3 #mS\n",
      "    APstim1.amp=1 #nA\n",
      "#Glutamate \n",
      "    Rel=h.STEP_REL(0.75,presyn)\n",
      "    Rel.amplitude=1\n",
      "    Rel.duration=3\n",
      "    Rel.release_time=120\n",
      "\n",
      "    for i in range(nofspines):\n",
      "        h.setpointer(Rel._ref_GLU,'C',AMPARsyn[i])\n",
      "        h.setpointer(Rel._ref_GLU,'C',NMDARsyn[i])\n",
      "    \n",
      "####Running the Simulation####\n",
      "    h.v_init=-85 #forced resting Vm for granule cells\n",
      "    h.init()\n",
      "\n",
      "    for l in range(nofsections):# dummy current source to compensate current caused by the forced Vm. \n",
      "        gc[l].ic_constant=-(gc[l].ina+gc[l].ik+gc[l].ica)\n",
      "    for l in range(nofspines):    \n",
      "        spineh[l].ic_constant=-(spineh[l].ina+spineh[l].ik+spineh[l].ica)\n",
      "\n",
      "    if h.cvode.active():\n",
      "        h.cvode.re_init()\n",
      "    else:\n",
      "        h.fcurrent()\n",
      "\n",
      "    h.tstop =350\n",
      "    h.run()\n",
      "    \n",
      "#### Vectors and conversion of units (um to nm)####\n",
      "    for i in range(nofspines):\n",
      "        v_spineh[i]=np.asarray(vrec_spineh[i])\n",
      "        ccaspineh_show[i]=1e6*np.asarray(ccaspineh[i])#converting to nM\n",
      "        icahvacc_show[i]=np.asarray(icahvacc[i]) #mA/cm2\n",
      "        icattype_show[i]=np.asarray(icattype[i]) #mA/cm2\n",
      "        ##to check the buffering capacity- not used in the figures\n",
      "        ccabuf1spineh_show[i]=1e6*np.asarray(ccabuf1spineh[i])\n",
      "        ccabuf2spineh_show[i]=1e6*np.asarray(ccabuf2spineh[i])\n",
      "        cbuf1spineh_show[i]=1e6*np.asarray(cbuf1spineh[i])\n",
      "        cbuf2spineh_show[i]=1e6*np.asarray(cbuf2spineh[i])\n",
      "        ##\n",
      "    for i in range(nofsections):\n",
      "        v_dend[i]=np.asarray(vrec_gc[i])\n",
      "        ccagc_show[i]=1e6*np.asarray(ccagc[i])#converting to nM\n",
      "\n",
      "####Mapping of Ca Concentration to df/f####\n",
      "    for i in range(nofspines): \n",
      "        y_dff_spineh[i]=spine_fit(ccaspineh_show[i])\n",
      "    \n",
      "    for i in range(nofsections):    \n",
      "        y_dff_gc[i]=spine_fit(ccagc_show[i])    \n",
      "       \n",
      "    pickle.dump(v_spineh[0],vmspine_obj)\n",
      "    pickle.dump(ccaspineh_show[0],caspine_0_obj)\n",
      "    pickle.dump(y_dff_spineh[0],dffspine_0_obj)\n",
      "    pickle.dump(ccagc_show[11],cagc_obj)\n",
      "    pickle.dump(y_dff_gc[11],dffgc_obj)\n",
      "    ##\n",
      "    pickle.dump(ccabuf1spineh_show[0],ccabuf1spineh_obj)\n",
      "    pickle.dump(ccabuf2spineh_show[0],ccabuf2spineh_obj)\n",
      "    pickle.dump(cbuf1spineh_show[0],cbuf1spineh_obj)\n",
      "    pickle.dump(cbuf2spineh_show[0],cbuf2spineh_obj)\n",
      "    ##\n",
      "    pickle.dump(icattype_show[0],ittype_obj)\n",
      "    pickle.dump(icahvacc_show[0],ihvacc_obj)\n",
      "  \n",
      "caspine_0_obj.close()\n",
      "dffspine_0_obj.close()\n",
      "vmspine_obj.close()\n",
      "cagc_obj.close()\n",
      "dffgc_obj.close()\n",
      "##\n",
      "ccabuf1spineh_obj.close()\n",
      "ccabuf2spineh_obj.close()\n",
      "cbuf1spineh_obj.close()\n",
      "cbuf2spineh_obj.close()\n",
      "##\n",
      "ittype_obj.close()\n",
      "ihvacc_obj.close()\n",
      "\n",
      "print \"run the next block\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "run the next block\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "####Retrieving the Simulation Data###\n",
      "####only the recorded variables which are used for SE are retrieving here.\n",
      "##cases for normalization\n",
      "t_c2=t_c=[-100,-80,-60,-40,-20,-12,-10,-8,-7,-6,-5,-4,-3,-2,0,2,3,4,5,6,7,8,10,12,20,40,60,80,100]\n",
      "t_c1=[0]\n",
      "\n",
      "y_dff_spineh_apalone=[np.array for i in range(len(t_c2))]\n",
      "y_dff_spineh_glualone=[np.array for i in range(len(t_c1))]\n",
      "cca_apalone=[np.array for i in range(len(t_c2))]\n",
      "cca_glualone=[np.array for i in range(len(t_c1))]\n",
      "vm_apalone=[np.array for i in range(len(t_c2))]\n",
      "vm_glualone=[np.array for i in range(len(t_c1))]\n",
      "\n",
      "dff_obj_glu=open(\"dffspine_0_record_glu\",\"r\")\n",
      "dff_obj_ap=open(\"dffspine_0_record_ap\",\"r\")\n",
      "cca_obj_glu=open(\"caspine_0_record_glu\",\"r\")\n",
      "cca_obj_ap=open(\"caspine_0_record_ap\",\"r\")\n",
      "vm_obj_ap=open(\"vmspine_0_record_ap\",\"r\")\n",
      "vm_obj_glu=open(\"vmspine_0_record_glu\",\"r\")\n",
      "\n",
      "for tc in range(len(t_c2)):\n",
      "    y_dff_spineh_apalone[tc]=pickle.load(dff_obj_ap)\n",
      "    cca_apalone[tc]=pickle.load(cca_obj_ap)\n",
      "    vm_apalone[tc]=pickle.load(vm_obj_ap)\n",
      " \n",
      "for tc in range(len(t_c1)):\n",
      "    y_dff_spineh_glualone[tc]=pickle.load(dff_obj_glu)\n",
      "    cca_glualone[tc]=pickle.load(cca_obj_glu)\n",
      "    vm_glualone[tc]=pickle.load(vm_obj_glu) \n",
      "\n",
      "dff_obj_glu.close()\n",
      "dff_obj_ap.close()\n",
      "cca_obj_ap.close()\n",
      "cca_obj_glu.close()\n",
      "vm_obj_glu.close()\n",
      "vm_obj_ap.close()\n",
      "\n",
      "##cases for coincidance\n",
      "ccaspineh_av_l=[np.array for i in range(len(t_c))]\n",
      "y_dff_spineh_av_l=[np.array for i in range(len(t_c))]\n",
      "vm_l=[np.array for i in range(len(t_c))]\n",
      "\n",
      "vm_obj=open(\"vmspine_0_record\",\"r\")\n",
      "caspine_0_obj=open(\"caspine_0_record\",\"r\")\n",
      "dffspine_0_obj=open(\"dffspine_0_record\",\"r\")\n",
      "\n",
      "for tc in range(len(t_c)):\n",
      "    ccaspineh_av_l[tc]=pickle.load(caspine_0_obj)\n",
      "    y_dff_spineh_av_l[tc]=pickle.load(dffspine_0_obj)\n",
      "    vm_l[tc]=pickle.load(vm_obj)\n",
      "\n",
      "caspine_0_obj.close()\n",
      "dffspine_0_obj.close()\n",
      "vm_obj.close()\n",
      "\n",
      "####Calculating SE####\n",
      "\n",
      "dff_add=[np.array for i in range(len(t_c))]\n",
      "cca_add=[np.array for i in range(len(t_c))]\n",
      "\n",
      "ext_ind_dff_add=[int for i in range(len(t_c))]\n",
      "ext_dff_add=[float for i in range(len(t_c))]\n",
      "\n",
      "ext_ind_cca_add=[int for i in range(len(t_c))]\n",
      "ext_cca_add=[float for i in range(len(t_c))]\n",
      "\n",
      "ext_ind_dff=[int for tc in range(len(t_c))]\n",
      "ext_dff=np.array([float for tc in range(len(t_c))])\n",
      "\n",
      "ext_ind_cca=[int for tc in range(len(t_c))]\n",
      "ext_cca=np.array([float for tc in range(len(t_c))])\n",
      "\n",
      "norm_cc=[float for i in range(len(t_c))]\n",
      "norm_dff=[float for i in range(len(t_c))]\n",
      "\n",
      "for tc in range(len(t_c)):\n",
      "\n",
      "    dff_add[tc]=y_dff_spineh_apalone[tc]+y_dff_spineh_glualone[0]\n",
      "    cca_add[tc]=cca_apalone[tc]+cca_glualone[0]-50\n",
      "\n",
      "for tc in range(len(t_c)):\n",
      "    ext_ind_dff_add[tc]=np.argmax(dff_add[tc])\n",
      "    ext_dff_add[tc]=dff_add[tc][ext_ind_dff_add[tc]]\n",
      "    \n",
      "    ext_ind_dff[tc]=np.argmax(y_dff_spineh_av_l[tc])\n",
      "    ext_dff[tc]=round(y_dff_spineh_av_l[tc][ext_ind_dff[tc]],3)\n",
      "    \n",
      "for tc in range(len(t_c)):\n",
      "    ext_ind_cca_add[tc]=np.argmax(cca_add[tc])\n",
      "    ext_cca_add[tc]=cca_add[tc][ext_ind_cca_add[tc]]\n",
      "    \n",
      "    ext_ind_cca[tc]=np.argmax(ccaspineh_av_l[tc])\n",
      "    ext_cca[tc]=round(ccaspineh_av_l[tc][ext_ind_cca[tc]],3)\n",
      "\n",
      "for tc in range(len(t_c)):    \n",
      "    norm_cc[tc]=round(ext_cca[tc]/ext_cca_add[tc],3)\n",
      "    norm_dff[tc]=round(ext_dff[tc]/ext_dff_add[tc],3)\n",
      "    \n",
      "####Sample Plot####\n",
      "plt.suptitle('norm',fontsize=14)\n",
      "plt.plot(t_c,norm_dff,'go')\n",
      "plt.grid()\n",
      "plt.show()\n",
      "\n",
      "print 'ext_dff=',ext_dff\n",
      "print 'norm_dff=',norm_dff"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "ext_dff= [64.533 65.295 65.938 65.756 61.584 56.311 54.302 51.753 50.13 48.035\n",
        " 38.531 35.885 36.972 38.327 41.211 43.005 44.1 46.573 50.403 54.694 58.321\n",
        " 61.096 65.234 68.378 76.271 81.193 80.246 78.278 76.106]\n",
        "norm_dff= [1.001, 0.998, 0.992, 0.972, 0.893, 0.81, 0.779, 0.741, 0.717, 0.687, 0.55, 0.512, 0.527, 0.546, 0.586, 0.611, 0.626, 0.66, 0.714, 0.775, 0.825, 0.864, 0.922, 0.966, 1.075, 1.146, 1.141, 1.126, 1.112]\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}