io_iv_curve.py

import json
import random
import time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import numba

@numba.jit(fastmath=True, cache=True, nopython=True)
def simulate(transient=10, seconds=1, delta=0.015, record_every=10):
    # Parameters
    g_int           =   0.13     # Cell internal conductance  -- now a parameter (0.13)
    p1              =   0.25     # Cell surface ratio soma/dendrite
    p2              =   0.15     # Cell surface ratio axon(hillock)/soma
    g_CaL           =   1.1      # Calcium T - (CaV 3.1) (0.7)
    g_h             =   0.12     # H current (HCN) (0.4996)
    g_K_Ca          =  35.0      # Potassium  (KCa v1.1 - BK) (35)
    g_ld            =   0.01532  # Leak dendrite (0.016)
    g_la            =   0.016    # Leak axon (0.016)
    g_ls            =   0.016    # Leak soma (0.016)
    S               =   1.0      # 1/C_m, cm^2/uF
    g_Na_s          = 150.0      # Sodium  - (Na v1.6 )
    g_Kdr_s         =   9.0      # Potassium - (K v4.3)
    g_K_s           =   5.0      # Potassium - (K v3.4)
    g_CaH           =   4.5      # High-threshold calcium -- Ca V2.1
    g_Na_a          = 240.0      # Sodium
    g_K_a           = 240.0      # Potassium (20)
    V_Na            =  55.0      # Sodium
    V_K             = -75.0      # Potassium
    V_Ca            = 120.0      # Low-threshold calcium channel
    V_h             = -43.0      # H current
    V_l             =  10.0      # Leak
    I_app           =   0.0

    g_CaL = 1

    # Soma state
    V_soma          = -60.0
    SOMA_k          =   0.7423159
    SOMA_l          =   0.0321349
    SOMA_h          =   0.3596066
    SOMA_n          =   0.2369847
    SOMA_x          =   0.1
    # Axon state
    V_axon          = -60.0
    AXON_Sodium_h   =   0.9
    AXON_Potassium_x=   0.2369847
    # Dend state
    V_dend          = -60.0
    DEND_Ca2Plus    =   3.715
    DEND_Calcium_r  =   0.0113
    DEND_Potassium_s=   0.0049291
    DEND_Hcurrent_q =   0.0337836

    def _update_soma(iv_trace, at):
        'Perform a single soma timestep update'
        nonlocal V_soma, SOMA_k, SOMA_l, SOMA_h, SOMA_n, SOMA_x
        # Leak current
        SOMA_I_leak        = g_ls * (V_soma - V_l)

        # Interaction Current
        I_ds        =  (g_int / p1)          * (V_soma - V_dend)
        I_as        =  (g_int / (1.0 - p2))  * (V_soma - V_axon)
        SOMA_I_interact =  I_ds + I_as

        # Channels
        # Low-threshold calcium
        SOMA_Ical   = g_CaL * SOMA_k * SOMA_k * SOMA_k * SOMA_l * (V_soma - V_Ca)

        SOMA_k_inf  = 1.0 / (1.0 + np.exp(-0.23809523809*V_soma - 14.5238))
        SOMA_l_inf  = 1.0 / (1.0 + np.exp( 0.11764705882*V_soma + 10.0))
        SOMA_tau_l  = (20.0 * np.exp(0.033333*V_soma + 5.333333) / (1.0  + np.exp(0.136986*V_soma + 11.506849))) + 35.0

        SOMA_dk_dt  = SOMA_k_inf - SOMA_k
        SOMA_dl_dt  = (SOMA_l_inf - SOMA_l) / SOMA_tau_l
        SOMA_k                  = delta * SOMA_dk_dt + SOMA_k
        SOMA_l                  = delta * SOMA_dl_dt + SOMA_l

        # Sodium   (watch out direct gate m)
        SOMA_m_inf  = 1.0 / (1.0 + np.exp(-0.1818181818*V_soma - 5.45454545))
        SOMA_Ina    = g_Na_s * SOMA_m_inf * SOMA_m_inf * SOMA_m_inf * SOMA_h * (V_soma - V_Na)
        SOMA_tau_h  =  3.0 * np.exp(0.0303030303*V_soma + 1.21212121212)
        SOMA_h_inf  = 1.0 / (1.0 + np.exp(0.1724137931*V_soma + 12.0689655))
        SOMA_dh_dt  = (SOMA_h_inf - SOMA_h) * SOMA_tau_h
        SOMA_h = SOMA_h + delta * SOMA_dh_dt

        # Potassium, slow component
        SOMA_Ikdr   = g_Kdr_s * SOMA_n * SOMA_n * SOMA_n * SOMA_n * (V_soma - V_K)
        SOMA_n_inf  = 1.0 / ( 1.0 + np.exp( -0.1*V_soma - 0.3))
        SOMA_tau_n  = 5.0 + (47.0 * np.exp(0.00111111*V_soma + 0.0555555555))
        SOMA_dn_dt  = SOMA_n_inf - SOMA_n / SOMA_tau_n
        SOMA_n = delta * SOMA_dn_dt + SOMA_n

        # Potassium, fast component
        SOMA_Ik      = g_K_s * SOMA_x**4 * (V_soma - V_K)
        SOMA_alpha_x = (0.13 * V_soma + 3.25) / (1.0 - np.exp(-0.1*V_soma - 2.5))
        SOMA_beta_x  = 1.69 * np.exp(-0.0125*V_soma -0.4375)
        SOMA_tau_x   = SOMA_alpha_x + SOMA_beta_x
        SOMA_x_inf   = SOMA_alpha_x / SOMA_tau_x

        SOMA_dx_dt   = (SOMA_x_inf - SOMA_x) * SOMA_tau_x
        SOMA_x  = delta * SOMA_dx_dt + SOMA_x

        if at >= 0:
            iv_trace[at, 0] = SOMA_Ik
            iv_trace[at, 1] = SOMA_Ikdr
            iv_trace[at, 2] = SOMA_Ina
            iv_trace[at, 3] = SOMA_Ical
            iv_trace[at, 4] = V_soma
        SOMA_I_Channels = SOMA_Ik + SOMA_Ikdr + SOMA_Ina + SOMA_Ical

        # Comp update
        SOMA_dv_dt = S * (-(SOMA_I_leak + SOMA_I_interact + SOMA_I_Channels))
        V_soma = V_soma + SOMA_dv_dt * delta

    def _update_axon(iv_trace, at):
        'Perform a single axon timestep update'
        nonlocal V_axon, AXON_Sodium_h, AXON_Potassium_x
        # Axon hillock components
        # Leak current
        AXON_I_leak    =  g_la * (V_axon - V_l)

        # Interaction Current
        I_sa           =  (g_int / p2) * (V_axon - V_soma)
        AXON_I_interact=   I_sa

        # Channelss
        # Sodium  (watch out direct gate !!!)
        AXON_m_inf     =  (1.0 / (1.0 + np.exp(-0.18181818*V_axon -5.45454545)))
        AXON_Ina       =  g_Na_a  * AXON_m_inf * AXON_m_inf * AXON_m_inf * AXON_Sodium_h * (V_axon - V_Na)
        AXON_h_inf     =  (1.0 / (1.0 + np.exp(0.1724137931*V_axon + 10.344827586)))
        AXON_tau_h     =  1.5 * np.exp(-0.0303030303*V_axon - 1.212121)
        AXON_dh_dt     =  ((AXON_h_inf - AXON_Sodium_h) /AXON_tau_h)
        AXON_Sodium_h              =  AXON_Sodium_h + delta * AXON_dh_dt

        # Potassium
        AXON_Ik        =  g_K_a * (AXON_Potassium_x)**4 * (V_axon - V_K)
        AXON_alpha_x   =  ((0.13*V_axon + 3.25) / (1.0 - np.exp(-0.1*V_axon - 2.5)))
        AXON_beta_x    =  1.69 * np.exp(-0.0125 * (V_axon + 35.0))
        AXON_x_inf     =  (AXON_alpha_x / (AXON_alpha_x + AXON_beta_x))
        AXON_tau_x     =  (1.0 / (AXON_alpha_x + AXON_beta_x))
        AXON_dx_dt     =  ((AXON_x_inf - AXON_Potassium_x) / AXON_tau_x)
        AXON_Potassium_x           =  delta * AXON_dx_dt + AXON_Potassium_x

        if at >= 0:
            iv_trace[at, 5] = AXON_Ina
            iv_trace[at, 6] = AXON_Ik
            iv_trace[at, 7] = V_axon
        AXON_I_Channels = AXON_Ina + AXON_Ik

        # comp update
        dv_dt  = S * (-(AXON_I_leak +  AXON_I_interact + AXON_I_Channels))
        V_axon = V_axon + dv_dt * delta

    def _update_dend(iv_trace, at):
        'Perform a single denrite timestep update'
        nonlocal V_dend, DEND_Ca2Plus, DEND_Calcium_r, DEND_Potassium_s, DEND_Hcurrent_q
        # Dendritic Components
        # Application current
        DEND_I_application = I_app

        # Leak current
        DEND_I_leak    =  g_ld * (V_dend - V_l)

        # Interaction Current
        DEND_I_interact        =  (g_int / (1.0 - p1)) * (V_dend - V_soma)

        # Channels
        # High-threshold calcium
        DEND_Icah        =  g_CaH * DEND_Calcium_r * DEND_Calcium_r * (V_dend - V_Ca)
        DEND_alpha_r     =  (1.7 / (1.0 + np.exp(-0.071942446*V_dend + 0.35971223021)))
        DEND_beta_r      =  (0.02*V_dend + 0.17) / (np.exp(0.2*V_dend + 1.7) - 1.0)
        DEND_tau_r       =  (DEND_alpha_r + DEND_beta_r)
        DEND_r_inf       =  (DEND_alpha_r / DEND_tau_r)
        DEND_dr_dt       =  (DEND_r_inf - DEND_Calcium_r) * DEND_tau_r * 0.2
        DEND_Calcium_r               =  delta * DEND_dr_dt + DEND_Calcium_r

        # Calcium dependent potassium
        DEND_Ikca        =  g_K_Ca * DEND_Potassium_s * (V_dend - V_K)
        DEND_alpha_s     =  (0.00002 * DEND_Ca2Plus) * (0.00002 * DEND_Ca2Plus < 0.01) + 0.01*((0.00002 * DEND_Ca2Plus)> 0.01)
        DEND_tau_s       =  DEND_alpha_s + 0.015
        DEND_s_inf       =  (DEND_alpha_s / DEND_tau_s)
        DEND_ds_dt       =  (DEND_s_inf - DEND_Potassium_s) * DEND_tau_s
        DEND_Potassium_s             =  delta * DEND_ds_dt + DEND_Potassium_s

        # calcium in general
        dCa_dt      =  -3.0 * DEND_Icah - 0.075 * DEND_Ca2Plus
        DEND_Ca2Plus            =  delta * dCa_dt + DEND_Ca2Plus

        # h current
        DEND_Ih     =  g_h * DEND_Hcurrent_q * (V_dend - V_h)
        q_inf       =  1.0 / (1.0 + np.exp(0.25*V_dend + 20.0))
        tau_q       =  np.exp(-0.086*V_dend - 14.6) + np.exp(0.070*V_dend - 1.87)
        dq_dt       =  (q_inf - DEND_Hcurrent_q) * tau_q
        DEND_Hcurrent_q         =  delta * dq_dt + DEND_Hcurrent_q

        if at >= 0:
            iv_trace[at, 8] = DEND_Icah
            iv_trace[at, 9] = DEND_Ikca
            iv_trace[at, 10] = DEND_Ih
            iv_trace[at, 11] = V_dend
        DEND_I_Channels = DEND_Icah + DEND_Ikca + DEND_Ih

        # comp update
        DEND_dv_dt  = S * (-(DEND_I_leak +  DEND_I_interact + DEND_I_application + DEND_I_Channels))
        V_dend = V_dend + DEND_dv_dt * delta

    #simulation
    nepochs = int(seconds*1000 / delta / record_every + .5)
    iv_trace = np.empty((nepochs, 12))
    nskip = int(1000 * transient / delta + 0.5)
    for i_skip in range(nskip):
        _update_soma(iv_trace, -1)
        _update_axon(iv_trace, -1)
        _update_dend(iv_trace, -1)
    for i_epoch in range(nepochs):
        for i_ts in range(record_every):
            _update_soma(iv_trace, i_epoch)
            _update_axon(iv_trace, i_epoch)
            _update_dend(iv_trace, i_epoch)
    return iv_trace


def main():
    iv_trace = simulate(10, 2)
    (SOMA_Ik, SOMA_Ikdr, SOMA_Ina, SOMA_Ical, V_soma,
     AXON_Ina, AXON_Ik, V_axon,
     DEND_Icah, DEND_Ikca, DEND_Ih, V_dend) = iv_trace.T

    plt.plot(V_soma)
    plt.plot(V_axon)
    plt.plot(V_dend)
    plt.show()


    # V_Na            =  55.0      # Sodium
    # V_K             = -75.0      # Potassium
    # V_Ca            = 120.0      # Low-threshold calcium channel
    # V_h             = -43.0      # H current

    fig, ax = plt.subplots(nrows=3, ncols=4, sharex=True, sharey=True)
    ax[0,0].set_title('Axon Ina')
    ax[0,0].plot(V_axon, AXON_Ina)
    ax[0,1].set_title('Axon Ik')
    ax[0,1].plot(V_axon, AXON_Ik)
    ax[0,2].axis('off')
    ax[0,3].axis('off')
    ax[1,0].set_title('Soma Ik')
    ax[1,0].plot(V_soma, SOMA_Ik)
    ax[1,1].set_title('Soma Ikdr')
    ax[1,1].plot(V_soma, SOMA_Ikdr)
    ax[1,2].set_title('Soma Ina')
    ax[1,2].plot(V_soma, SOMA_Ina)
    ax[1,3].set_title('Soma Ical')
    ax[1,3].plot(V_soma, SOMA_Ical)
    ax[2,0].set_title('Dend Icah')
    ax[2,0].plot(V_dend, DEND_Icah)
    ax[2,1].set_title('Dend Ikca')
    ax[2,1].plot(V_dend, DEND_Ikca)
    ax[2,2].set_title('Dend Ih')
    ax[2,2].plot(V_dend, DEND_Ih)
    ax[2,3].axis('off')
    #
    ax[2,0].set_xlabel('V (mV)')
    ax[2,1].set_xlabel('V (mV)')
    ax[2,2].set_xlabel('V (mV)')
    ax[1,3].set_xlabel('V (mV)')
    #
    ax[0,0].set_ylabel('I')
    ax[1,0].set_ylabel('I')
    ax[2,0].set_ylabel('I')
    #
    plt.tight_layout()
    plt.show()


if __name__ == '__main__':
    main()