{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# import packages and define stimulus protocols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# using MAT\n",
    "\n",
    "using MultivariateStats\n",
    "using DifferentialEquations\n",
    "using Plots\n",
    "gr(fmt=\"png\", size=(1600, 600))\n",
    "\n",
    "include(\"../src/HHModel.jl\")\n",
    "\n",
    "# stimulus protocol\n",
    "current_step = (t, param) -> begin\n",
    "    (param.start < t)&(t < param.dur + param.start) ? param.step + param.holding + param.noise * randn() : param.holding + param.holding + param.noise * randn()\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## setting up biophysical parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type one in matlab - sustained\n",
    "\n",
    "# biophysical model setup ==> usu. for maximum conductance\n",
    "htk = HHModel.high_voltage_gated_potassium(2.8, phi=0.85)\n",
    "# ltk = HHModel.low_voltage_gated_potassium(0.65, subtype=:kv1)\n",
    "# ik = HHModel.hh_potassium(10.0)\n",
    "ina = HHModel.hh_sodium(20.0)\n",
    "ih = HHModel.ihcurrent(0.91)\n",
    "il = HHModel.leakage(0.03)\n",
    "\n",
    "_model = [htk, ina, ih, il]\n",
    "_model_sim_cc = HHModel.simpleConductanceModel(_model, current_step, C=0.9); # current clamp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# type two in matlab -- transient\n",
    "\n",
    "# biophysical model setup ==> usu. for maximum conductance\n",
    "htk = HHModel.high_voltage_gated_potassium(2.8, phi=0.85)\n",
    "ltk = HHModel.low_voltage_gated_potassium(1.1, subtype=:kv1)\n",
    "# ik = HHModel.hh_potassium(10.0)\n",
    "ina = HHModel.hh_sodium(13.0)\n",
    "ih = HHModel.ihcurrent(0.43)\n",
    "il = HHModel.leakage(0.03)\n",
    "\n",
    "_model = [htk, ltk, ina, ih, il]\n",
    "_model_sim_cc = HHModel.simpleConductanceModel(_model, current_step, C=0.9); # current clamp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# custom parameters\n",
    "# biophysical model setup ==> usu. for maximum conductance\n",
    "htk = HHModel.high_voltage_gated_potassium(1.8)\n",
    "ltk = HHModel.low_voltage_gated_potassium(0.1, subtype=:kv1)\n",
    "# ik = HHModel.hh_potassium(10.0)\n",
    "ina = HHModel.hh_sodium(13.0)\n",
    "ih = HHModel.ihcurrent(0.3)\n",
    "il = HHModel.leakage(0.2)\n",
    "\n",
    "_model = [htk, ltk, ina, ih, il]\n",
    "# _model = [htk, ina, ih, il]\n",
    "_model_sim_cc = HHModel.simpleConductanceModel(_model, current_step, C=0.9); # current clamp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## running the simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# running the simulation\n",
    "_p = (E=(sodium=82.0, potassium=-81.0, ih=-46.0, leak=-65.0, syn=3.0), \n",
    "      stim=(start=500, step=5.0, holding=0, dur=1500, noise=0),)\n",
    "\n",
    "# ih.h.Vhalf = -60.0\n",
    "# HHModel.update!(ih)\n",
    "\n",
    "tspan = (0.0, 2500.0)\n",
    "v0 = -60.9766\n",
    "v0 = -72.7030;\n",
    "u0 = HHModel.setup_init(_model, v0)\n",
    "\n",
    "prob = ODEProblem(_model_sim_cc, u0, tspan, _p)\n",
    "@time sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)\n",
    "\n",
    "# preview voltage trace\n",
    "plot(sol, vars=(1), legend=nothing)\n",
    "# savefig(\"demo2.svg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# preview current\n",
    "_trange = 400:0.1:1000\n",
    "_trace = HHModel.current_decompose(sol, _model, _trange, _p)\n",
    "plot(xlim=(_trange[1], _trange[end]), legend=:right)\n",
    "for (key, val) in _trace\n",
    "    if key == \"voltage\"\n",
    "        continue\n",
    "    end\n",
    "    plot!(_trange, val, label=key, linewidth=3)\n",
    "end\n",
    "plot!()\n",
    "# savefig(\"demo2_current.svg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# export\n",
    "_export = HHModel.current_decompose(sol, _model, 0.01:0.01:2500, _p)\n",
    "\n",
    "# _export = Dict(\"param\" => hcat(sol(0.01:0.01:2500).u...)[1, :])\n",
    "MAT.matwrite(\"julia_sim_transient_currents.mat\", _export)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# preview dynamics\n",
    "\n",
    "# dynamics between membrane potential and a particular kinetic variable\n",
    "plot(sol, vars=(1,2), size=(500, 500), legend=nothing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_down_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PCA projection\n",
    "_var_array = hcat(sol(500:0.1:2000).u...)[2:end-1, :]\n",
    "\n",
    "_down = fit(PCA, _var_array, );\n",
    "\n",
    "_down_array = transform(_down, _var_array);\n",
    "\n",
    "plot(_down_array[1,:], _down_array[2, :], size=(500, 500), legend=nothing)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sweep protocol\n",
    "\n",
    "_current_steps = -5:0.5:0\n",
    "\n",
    "plot(legend=nothing, ylim=(-150, 80), xlim=(190, 250))\n",
    "\n",
    "@time for _i_step in _current_steps\n",
    "    _p = (E=(sodium=81.27, potassium=-80.78, ih=-26, leak=-65, syn=3), \n",
    "      stim=(start=50, step=_i_step, dur=500, holding=0, noise=0),)\n",
    "    prob = ODEProblem(_model_sim_cc, u0, tspan, _p)\n",
    "    sol = solve(prob, Tsit5(), dt=0.1)\n",
    "\n",
    "    plot!(sol, vars=(1), color=:black, alpha=0.3)\n",
    "end\n",
    "\n",
    "plot!()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# voltage clamp\n",
    "voltage_step = (t, param) -> begin\n",
    "    (param.start < t)&(t < param.dur + param.start) ? param.step + param.noise * (rand()-0.5) : -60 + param.noise * (rand()-0.5)\n",
    "end\n",
    "\n",
    "_model_sim = HHModel.simpleVoltageClamp(_model, voltage_step, C=0.9)\n",
    "\n",
    "# running the simulation\n",
    "_p = (E=(sodium=81.27, potassium=-80.78, ih=-26.0, leak=-65.0, syn=3.0), \n",
    "      stim=(start=200, step=30, dur=1500, noise=0),)\n",
    "\n",
    "tspan = (0.0, 2000.0)\n",
    "v0 = -60.0\n",
    "u0 = HHModel.setup_init(_model, v0)\n",
    "\n",
    "prob = ODEProblem(_model_sim, u0, tspan, _p)\n",
    "@time sol = solve(prob, Tsit5(), dt=0.01)\n",
    "\n",
    "# preview current\n",
    "plot(sol, vars=(length(u0)), legend=nothing, xlim=(180,250))\n",
    "# savefig(\"demo2.svg\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# preview current\n",
    "_trange = 0:0.1:2000\n",
    "_trace = HHModel.current_decompose(sol, _model, _trange, _p)\n",
    "plot(xlim=(_trange[1], _trange[end]))\n",
    "for (key, val) in _trace\n",
    "    if key == \"voltage\"\n",
    "        continue\n",
    "    end\n",
    "    plot!(_trange, val, label=key)\n",
    "end\n",
    "plot!(xlim=(180, 1000))\n",
    "# savefig(\"demo2_current.svg\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.2.0",
   "language": "julia",
   "name": "julia-1.2"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}