MakeFig1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fig 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using PyPlot\n",
    "using Statistics\n",
    "using Distributions   \n",
    "using Random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Functions  \n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "function EIFshotEI0I(T,dt,vth,vre,vlb,tau,vT,dT,ee,ei,se,si)\n",
    "\n",
    "    n=Int(T/dt)\n",
    "    t=dt*collect(1:n)\n",
    "    \n",
    "    # initialize and generate random numbers\n",
    "    c=0\n",
    "    v=zeros(n)\n",
    "    tp=zeros(Int,n)\n",
    "    \n",
    "    for k=1:n-1\n",
    "        v[k+1]=v[k] + (dt/tau)*(dT*exp((v[k]-vT)/dT)-v[k]) + (ee-0)*se[k]+(ei-0)*si[k];\n",
    "        v[k+1]>vth ? (c=c+1; v[k]=vth; v[k+1]=vre; tp[k]=1) : ()\n",
    "        v[k+1]<vlb ?  v[k+1]=vlb : ()\n",
    "    end\n",
    "    \n",
    "    r=c/T\n",
    "    tp=dt*findall(tp.==1) # the spike time  \n",
    "    \n",
    "    return t,v,r,tp\n",
    "        \n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "function EIFshotEI0g(T,dt,vth,vre,vlb,tau,vT,dT,ee,ei,se,si)\n",
    "\n",
    "    n=Int(T/dt)\n",
    "    t=dt*collect(1:n)\n",
    "    \n",
    "    # initialize and generate random numbers\n",
    "    c=0\n",
    "    v=zeros(n)\n",
    "    tp=zeros(Int,n)\n",
    "    \n",
    "    for k=1:n-1\n",
    "        v[k+1]=v[k] + (dt/tau)*(dT*exp((v[k]-vT)/dT)-v[k]) + (ee-v[k])*se[k] + (ei-v[k])*si[k] \n",
    "        v[k+1]>vth ? (c=c+1; v[k]=vth; v[k+1]=vre; tp[k]=1) : ()\n",
    "        v[k+1]<vlb ?  v[k+1]=vlb : ()\n",
    "    end\n",
    "    \n",
    "    r=c/T\n",
    "    tp=dt*findall(tp.==1) # the spike time\n",
    "    \n",
    "    return t,v,r,tp\n",
    "\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Main code   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Neuron params:\ttau=20ms, ee=60mV, ei=-10mV vT=10mV, dT=1mV\n",
      "Thresh params:\tvthEIF=20mV vthLIF=10\n",
      "LowerB params:\tvlbCond=-9.9mV vlbCurr=-15.0\n",
      "Syn params:\tae,ai=1.5, -0.75mV and sig=5mV\n",
      "Num params:\tThIn dv=0.0005mV & sims dt=0.01ms\n"
     ]
    }
   ],
   "source": [
    "include(\"CommonCode.jl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "###########################################\n",
    "# Panel A - forcing term\n",
    "###########################################\n",
    "\n",
    "# get the LIF forcing t erm\n",
    "vvLIF=collect(-10:0.1:vthLIF)\n",
    "FLIF=-vvLIF/tau\n",
    "FLIF[end]=10\n",
    "\n",
    "# get the EIF forcing terms\n",
    "vv=collect(-10:0.1:20)\n",
    "FEIF=(-vv .+ dT*exp.((vv.-vT)/dT))/tau;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "###########################################\n",
    "# Panel B - distributions\n",
    "###########################################\n",
    "\n",
    "# parameters for current-based input\n",
    "a=collect(-10:0.1:10)\n",
    "\n",
    "Ae=(a.>0).*exp.(-a/ae)/ae\n",
    "Ai=-(a.<0).*exp.(-a/ai)/ai\n",
    "\n",
    "# parameters for conductance-based input\n",
    "\n",
    "b=collect(0:0.01:1)\n",
    "\n",
    "be=ae/ee\n",
    "bi=ai/ei\n",
    "betae=1/be-1\n",
    "betai=1/bi-1\n",
    "Be=betae*(1 .-b).^(betae-1)\n",
    "Bi=betai*(1 .-b).^(betai-1)\n",
    "\n",
    "v1=-5\n",
    "aeg1=v1 .+b*(ee-v1); Aeg1=Be/(ee-v1)\n",
    "aig1=v1 .+b*(ei-v1); Aig1=-Bi/(ei-v1)\n",
    "\n",
    "v2=0;\n",
    "aeg2=v2 .+b*(ee-v2); Aeg2=Be/(ee-v2)\n",
    "aig2=v2 .+b*(ei-v2); Aig2=-Bi/(ei-v2)\n",
    "\n",
    "v3=5;\n",
    "aeg3=v3 .+b*(ee-v3); Aeg3=Be/(ee-v3)\n",
    "aig3=v3 .+b*(ei-v3); Aig3=-Bi/(ei-v3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "###########################################\n",
    "# panel C - timeseries examples\n",
    "###########################################\n",
    "\n",
    "Random.seed!(10);\n",
    "\n",
    "T=500\n",
    "nt=Int(T/dt);\n",
    "\n",
    "Re,Ri=0.2,0.02\n",
    "\n",
    "# Specific EIF PARAMS\n",
    "\n",
    "be=ae/ee;\n",
    "bi=ai/ei;\n",
    "betae=1/be-1;\n",
    "betai=1/bi-1;\n",
    "\n",
    "se=rand(Bernoulli(Re*dt),nt).*rand(Beta(1,betae),nt) \n",
    "si=rand(Bernoulli(Ri*dt),nt).*rand(Beta(1,betai),nt) \n",
    "\n",
    "# LIF sims for current and conductance\n",
    "vthLIF=10;\n",
    "t,vLIFI,rLIFI,tpLIFI=EIFshotEI0I(T,dt,vthLIF,vre,vlbCurr,tau,vT,0.0,ee,ei,se,si)\n",
    "t,vLIFg,rLIFg,tpLIFg=EIFshotEI0g(T,dt,vthLIF,vre,vlbCond,tau,vT,0.0,ee,ei,se,si)\n",
    "\n",
    "# EIF sims for current and conductance\n",
    "vthEIF=50\n",
    "t,vEIFI,rEIFI,tpEIFI=EIFshotEI0I(T,dt,vthEIF,vre,vlbCurr,tau,vT,dT,ee,ei,se,si)\n",
    "t,vEIFg,rEIFg,tpEIFg=EIFshotEI0g(T,dt,vthEIF,vre,vlbCond,tau,vT,dT,ee,ei,se,si);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "##################################################################\n",
    "# Create Figure 1\n",
    "##################################################################\n",
    "\n",
    "pygui(true); fig=figure(figsize=1.6.*(3.375,3))    \n",
    "\n",
    "##################################################################\n",
    "# Positions\n",
    "##################################################################\n",
    "bx,tx=0.1,0.95 # top and bottom x edges\n",
    "by,ty=0.1,0.925 # top and bottom y edges\n",
    "gx,gy=0.12,0.12 # gap sizes\n",
    "\n",
    "panelfs=10\n",
    "labelfs=10\n",
    "tickfs=10\n",
    "lw=1.0\n",
    "ms=3\n",
    "\n",
    "EIFCol=\"black\"\n",
    "LIFCol=\"black\"\n",
    "EIFCondCol=\"black\"\n",
    "EIFCurrCol=\"darkgrey\"\n",
    "LIFCondCol=\"black\"\n",
    "LIFCurrCol=\"darkgrey\"\n",
    "\n",
    "##################################################################\n",
    "##################################################################\n",
    "# PANEL 1A\n",
    "##################################################################\n",
    "##################################################################\n",
    "xA,yA=0.1,0.525\n",
    "wA,hA=0.375,0.4\n",
    "\n",
    "axA=plt.axes([bx,yA,wA,hA]);\n",
    "fig.text(xA-0.0925, 0.95, \"(a)\", fontsize=panelfs)\n",
    "plot(vvLIF,FLIF,\"-\",color=LIFCol,linewidth=lw)\n",
    "plot(vv,FEIF,\"-\",color=EIFCol,linewidth=lw)\n",
    "plot([-10,20],[0,0],\"k:\",linewidth=lw)\n",
    "xlabel(\"v (mv)\",fontsize=labelfs); \n",
    "ylabel(\"f(v) (mV/ms)\",fontsize=labelfs);\n",
    "xticks(fontsize=tickfs);yticks(fontsize=tickfs)\n",
    "axis([-10,15,-1,2]);\n",
    "axA.set_yticks([-1,0,1,2]);\n",
    "axA.set_yticklabels([\"-1\",\"0\",\"1\",\"2\"],fontsize=tickfs)\n",
    "axA.text(7, 1.5, \"LIF\", fontsize=panelfs)\n",
    "axA.text(11, 1.5, \"EIF\", fontsize=panelfs)\n",
    "\n",
    "##################################################################\n",
    "##################################################################\n",
    "# PANEL 1B\n",
    "##################################################################\n",
    "##################################################################\n",
    "xB,yB1,yB2=0.6,0.525,0.8\n",
    "wB,hB=0.35,0.125\n",
    "\n",
    "fig.text(xB-0.115, 0.95, \"(b)\", fontsize=panelfs)\n",
    "\n",
    "axB2=plt.axes([xB,yB2,wB,hB]);\n",
    "EPSPcolor=\"green\"\n",
    "IPSPcolor=\"red\"\n",
    "fill_between(aeg1, Aeg1, 0, color=EPSPcolor, alpha=0.3)\n",
    "fill_between(aig1, Aig1, 0, color=IPSPcolor, alpha=0.3)\n",
    "fill_between(aeg2, Aeg2, 0, color=EPSPcolor, alpha=0.3)\n",
    "fill_between(aig2, Aig2, 0, color=IPSPcolor, alpha=0.3)\n",
    "fill_between(aeg3, Aeg3, 0, color=EPSPcolor, alpha=0.3)\n",
    "fill_between(aig3, Aig3, 0, color=IPSPcolor, alpha=0.3)\n",
    "axis([-10,10,0,2.5])\n",
    "ylabel(\"dist (cond)\",fontsize=labelfs); \n",
    "\n",
    "fig.text(0.725, 0.715, \"from v=0\", fontsize=tickfs)\n",
    "fig.text(0.725, 0.68, L\"a_\\mathrm{e}=1.5\\mathrm{mV}\", fontsize=tickfs)\n",
    "fig.text(0.725, 0.645, L\"a_\\mathrm{i}=-0.75\\mathrm{mV}\", fontsize=tickfs)\n",
    "\n",
    "axB2.set_yticks([0,3]);\n",
    "axB2.set_yticklabels([\"0\",\"3\"],fontsize=tickfs)\n",
    "#yticks(rotation = 90,fontsize=tickfs)\n",
    "\n",
    "axB1=plt.axes([xB,yB1,wB,hB]);\n",
    "fill_between(a .+v1, Ae, 0, color=EPSPcolor, alpha=0.3)\n",
    "fill_between(a .+v1, Ai, 0, color=IPSPcolor, alpha=0.3)\n",
    "fill_between(a .+v2, Ae, 0, color=EPSPcolor, alpha=0.3)\n",
    "fill_between(a .+v2, Ai, 0, color=IPSPcolor, alpha=0.3)\n",
    "fill_between(a .+v3, Ae, 0, color=EPSPcolor, alpha=0.3)\n",
    "fill_between(a .+v3, Ai, 0, color=IPSPcolor, alpha=0.3)\n",
    "axis([-10,10,0,2.5]);\n",
    "ylabel(\"dist (curr)\",fontsize=labelfs);\n",
    "xlabel(\"Post-synaptic potentials (mV)\",fontsize=labelfs);\n",
    "#xticks(fontsize=tickfs);yticks(fontsize=tickfs)\n",
    "axB1.set_yticks([0,3]);\n",
    "axB1.set_yticklabels([\"0\",\"3\"],fontsize=tickfs)\n",
    "\n",
    "##################################################################\n",
    "##################################################################\n",
    "# PANEL 1C\n",
    "##################################################################\n",
    "##################################################################\n",
    "xC=0.1\n",
    "yC1,yC2=0.275,0.1\n",
    "wC,hC=0.85,0.125\n",
    "\n",
    "LIFspike=14\n",
    "EIFspike=16\n",
    "\n",
    "fig.text(xC-0.0925, 0.425, \"(c)\", fontsize=panelfs)\n",
    "\n",
    "axC1=plt.axes([bx,yC1,wC,hC]);\n",
    "plot(t,vLIFI,color=LIFCurrCol,linewidth=lw)\n",
    "plot(t,vLIFg,color=LIFCondCol,linewidth=lw)\n",
    "plot(tpLIFI,LIFspike .+ zeros(length(tpLIFI)),\"*\",color=LIFCurrCol)\n",
    "plot(tpLIFg,LIFspike .+ zeros(length(tpLIFg)),\"*\",color=LIFCondCol)\n",
    "ylabel(\"v (mV)\",fontsize=labelfs)\n",
    "xticks(fontsize=tickfs);yticks(fontsize=tickfs)\n",
    "axis([0,500,0,17])\n",
    "axC1.text(70, 10, \"LIF\", fontsize=panelfs)\n",
    "\n",
    "axC2=plt.axes([bx,yC2,wC,hC]);\n",
    "plot(t,vEIFI,color=EIFCurrCol,linewidth=lw)\n",
    "plot(t,vEIFg,color=EIFCondCol,linewidth=lw)\n",
    "plot(tpEIFI,EIFspike .+zeros(length(tpEIFI)),\"*\",color=EIFCurrCol)\n",
    "plot(tpEIFg,EIFspike .+zeros(length(tpEIFg)),\"*\",color=EIFCondCol)\n",
    "ylabel(\"v (mV)\",fontsize=labelfs); xlabel(\"Time (ms)\",fontsize=labelfs)\n",
    "xticks(fontsize=tickfs);yticks(fontsize=tickfs)\n",
    "axis([0,500,0,17])\n",
    "axC2.text(70, 10, \"EIF\", fontsize=panelfs)\n",
    "\n",
    "##################################################################\n",
    "##################################################################\n",
    "# TIDY UP ALL FIGS\n",
    "##################################################################\n",
    "##################################################################\n",
    "for ax in [axA,axB1,axB2,axC1,axC2]\n",
    "    ax.spines[\"top\"].set_visible(false)\n",
    "    ax.spines[\"right\"].set_visible(false)\n",
    "    sca(ax); yticks(fontsize=tickfs); xticks(fontsize=tickfs); \n",
    "end\n",
    "\n",
    "savefig(\"Fig1.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.4.1",
   "language": "julia",
   "name": "julia-1.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.4.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}