simulator_tests.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95b93c7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# STEPS test(https://github.com/CNS-OIST/STEPS_Example/blob/master/python_scripts/API_2/ip3/ip3.py)\n",
    "import steps.interface\n",
    "\n",
    "from steps.model import *\n",
    "from steps.geom import *\n",
    "from steps.sim import *\n",
    "from steps.saving import *\n",
    "from steps.rng import *\n",
    "\n",
    "import sys\n",
    "import pylab\n",
    "import numpy\n",
    "\n",
    "r = ReactionManager()\n",
    "\n",
    "#########################\n",
    "# Model setup\n",
    "#########################\n",
    "mdl = Model()\n",
    "with mdl:\n",
    "    # Species \n",
    "    Ca, IP3, R, RIP3, Ropen, RCa, R2Ca, R3Ca, R4Ca = Species.Create()\n",
    "    # Surface system\n",
    "    surfsys = SurfaceSystem.Create()\n",
    "    # Reactions\n",
    "    with surfsys:\n",
    "        # IP3 binding\n",
    "        (IP3.o + R.s <r[1]> RIP3.s) + Ca.o <r[2]> Ropen.s\n",
    "        r[1].K = 1000e6, 25800\n",
    "        r[2].K = 8000e6, 2000\n",
    "\n",
    "        # Ca binding\n",
    "        (((R.s + Ca.o <r[7]> RCa.s) + Ca.o <r[8]> R2Ca.s) + Ca.o <r[3]> R3Ca.s) + Ca.o <r[4]> R4Ca.s\n",
    "        r[7].K = 8.889e6, 5\n",
    "        r[8].K = 20e6, 10\n",
    "        r[3].K = 40e6, 15\n",
    "        r[4].K = 60e6, 20\n",
    "\n",
    "        # Ca ions passing through open IP3R channel\n",
    "        Ca.i + Ropen.s >r[1]> Ropen.s + Ca.o\n",
    "        r[1].K = 2e8\n",
    "\n",
    "#########################\n",
    "# Geom setup\n",
    "#########################\n",
    "geom = Geometry()\n",
    "with geom:\n",
    "    # Create the cytosol and Endoplasmic Reticulum compartments\n",
    "    cyt, ER = Compartment.Create()\n",
    "    cyt.Vol = 1.6572e-19\n",
    "    ER.Vol = 1.968e-20\n",
    "\n",
    "    # ER is the 'inner' compartment, cyt is the 'outer' compartment\n",
    "    memb = Patch.Create(ER, cyt, surfsys)\n",
    "    memb.Area = 0.4143e-12\n",
    "\n",
    "    print('Inner compartment to memb is', memb.innerComp.name)\n",
    "    print('Outer compartment to memb is', memb.outerComp.name)\n",
    "\n",
    "#########################\n",
    "# RNG setup\n",
    "#########################\n",
    "rng = RNG('mt19937', 512, 7233)\n",
    "\n",
    "#########################\n",
    "# Solver setup\n",
    "#########################\n",
    "sim = Simulation('Wmdirect', mdl, geom, rng)\n",
    "        \n",
    "NITER = 100\n",
    "ENDT = 0.201\n",
    "\n",
    "rs = ResultSelector(sim)\n",
    "\n",
    "RopenDat = rs.memb.Ropen.Count\n",
    "\n",
    "CaConc = rs.cyt.Ca.Conc\n",
    "\n",
    "sim.toSave(RopenDat, CaConc, dt=0.001)\n",
    "\n",
    "#########################\n",
    "# Run simulation\n",
    "#########################\n",
    "\n",
    "for i in range (0, NITER):\n",
    "    sim.newRun()\n",
    "\n",
    "    sim.cyt.Ca.Conc = 3.30657e-8\n",
    "    sim.cyt.IP3.Count = 6\n",
    "    sim.ER.Ca.Conc = 150e-6\n",
    "    sim.ER.Ca.Clamped = True\n",
    "    sim.memb.R.Count = 160\n",
    "\n",
    "    sim.run(ENDT)\n",
    "\n",
    "    pylab.plot(RopenDat.time[-1], RopenDat.data[-1], color = 'blue', linewidth = 0.1)\n",
    "\n",
    "datTime = RopenDat.time[0]\n",
    "\n",
    "res_mean = numpy.mean(RopenDat.data, 0)\n",
    "res_std = numpy.std(RopenDat.data, 0)\n",
    "res_std1 = res_mean + res_std\n",
    "res_std2 = res_mean - res_std\n",
    "\n",
    "pylab.plot(datTime, res_mean, color = 'black', linewidth = 2.0, label = 'mean')\n",
    "pylab.plot(datTime, res_std1, color = 'gray', linewidth = 1.0, label='std')\n",
    "pylab.plot(datTime, res_std2,color = 'gray', linewidth = 1.0)\n",
    "\n",
    "pylab.xlabel('Time (sec)')\n",
    "pylab.ylabel('# IP3 receptors in open state')\n",
    "pylab.title('IP3 receptor model: %d iterations with Wmdirect'%NITER)\n",
    "pylab.ylim(0)\n",
    "pylab.legend()\n",
    "\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53693961",
   "metadata": {},
   "outputs": [],
   "source": [
    "#NEURON test (OCNC 2019 NEURON Tutorial 1)\n",
    "from neuron import h\n",
    "h.load_file('stdrun.hoc')\n",
    "import matplotlib.pyplot as plt\n",
    "soma = h.Section(name='soma')\n",
    "soma.insert('pas')\n",
    "asyn = h.AlphaSynapse(soma(0.5))\n",
    "asyn.onset = 20.0\n",
    "asyn.gmax = 1.0\n",
    "asyn.tau = 0.1\n",
    "asyn.e = 0.0\n",
    "v_vec = h.Vector()              \n",
    "t_vec = h.Vector()              \n",
    "v_vec.record(soma(0.5)._ref_v)\n",
    "t_vec.record(h._ref_t)\n",
    "h.tstop = 40.0\n",
    "h.run()\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.plot(t_vec, v_vec)\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('Soma membrane potential (mV)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3154cd24",
   "metadata": {},
   "outputs": [],
   "source": [
    "#NEST test (https://nest-simulator.readthedocs.io/en/stable/auto_examples/one_neuron.html)\n",
    "import nest\n",
    "import nest.voltage_trace\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "nest.set_verbosity(\"M_WARNING\")\n",
    "nest.ResetKernel()\n",
    "neuron = nest.Create(\"iaf_psc_alpha\")\n",
    "voltmeter = nest.Create(\"voltmeter\")\n",
    "neuron.I_e = 376.0\n",
    "nest.Connect(voltmeter, neuron)\n",
    "nest.Simulate(1000.0)\n",
    "nest.voltage_trace.from_device(voltmeter)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b95a1943",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "237.4 ms (23%) simulated in 3m 29s, estimated 11m 13s remaining.\n"
     ]
    }
   ],
   "source": [
    "#Brian2 test (https://brian2.readthedocs.io/en/stable/examples/COBAHH.html)\n",
    "from brian2 import *\n",
    "\n",
    "# Parameters\n",
    "area = 20000*umetre**2\n",
    "Cm = (1*ufarad*cm**-2) * area\n",
    "gl = (5e-5*siemens*cm**-2) * area\n",
    "\n",
    "El = -60*mV\n",
    "EK = -90*mV\n",
    "ENa = 50*mV\n",
    "g_na = (100*msiemens*cm**-2) * area\n",
    "g_kd = (30*msiemens*cm**-2) * area\n",
    "VT = -63*mV\n",
    "# Time constants\n",
    "taue = 5*ms\n",
    "taui = 10*ms\n",
    "# Reversal potentials\n",
    "Ee = 0*mV\n",
    "Ei = -80*mV\n",
    "we = 6*nS  # excitatory synaptic weight\n",
    "wi = 67*nS  # inhibitory synaptic weight\n",
    "\n",
    "# The model\n",
    "eqs = Equations('''\n",
    "dv/dt = (gl*(El-v)+ge*(Ee-v)+gi*(Ei-v)-\n",
    "         g_na*(m*m*m)*h*(v-ENa)-\n",
    "         g_kd*(n*n*n*n)*(v-EK))/Cm : volt\n",
    "dm/dt = alpha_m*(1-m)-beta_m*m : 1\n",
    "dn/dt = alpha_n*(1-n)-beta_n*n : 1\n",
    "dh/dt = alpha_h*(1-h)-beta_h*h : 1\n",
    "dge/dt = -ge*(1./taue) : siemens\n",
    "dgi/dt = -gi*(1./taui) : siemens\n",
    "alpha_m = 0.32*(mV**-1)*4*mV/exprel((13*mV-v+VT)/(4*mV))/ms : Hz\n",
    "beta_m = 0.28*(mV**-1)*5*mV/exprel((v-VT-40*mV)/(5*mV))/ms : Hz\n",
    "alpha_h = 0.128*exp((17*mV-v+VT)/(18*mV))/ms : Hz\n",
    "beta_h = 4./(1+exp((40*mV-v+VT)/(5*mV)))/ms : Hz\n",
    "alpha_n = 0.032*(mV**-1)*5*mV/exprel((15*mV-v+VT)/(5*mV))/ms : Hz\n",
    "beta_n = .5*exp((10*mV-v+VT)/(40*mV))/ms : Hz\n",
    "''')\n",
    "\n",
    "P = NeuronGroup(4000, model=eqs, threshold='v>-20*mV', refractory=3*ms,\n",
    "                method='exponential_euler')\n",
    "Pe = P[:3200]\n",
    "Pi = P[3200:]\n",
    "Ce = Synapses(Pe, P, on_pre='ge+=we')\n",
    "Ci = Synapses(Pi, P, on_pre='gi+=wi')\n",
    "Ce.connect(p=0.02)\n",
    "Ci.connect(p=0.02)\n",
    "\n",
    "# Initialization\n",
    "P.v = 'El + (randn() * 5 - 5)*mV'\n",
    "P.ge = '(randn() * 1.5 + 4) * 10.*nS'\n",
    "P.gi = '(randn() * 12 + 20) * 10.*nS'\n",
    "\n",
    "# Record a few traces\n",
    "trace = StateMonitor(P, 'v', record=[1, 10, 100])\n",
    "run(1 * second, report='text')\n",
    "plot(trace.t/ms, trace[1].v/mV)\n",
    "plot(trace.t/ms, trace[10].v/mV)\n",
    "plot(trace.t/ms, trace[100].v/mV)\n",
    "xlabel('t (ms)')\n",
    "ylabel('v (mV)')\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c82ca82c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}