Activity-mediated_accumulation_potassium_induces_switch_firing_pattern/cfg/m_MTM_W_sPNaS_sICD.py

from scipy import integrate
from scipy.integrate import odeint
import numpy as np
from copy import copy

Pars_MTM_W_sPNaS_sICD={
    'I_app':1, #uA/cm2
    'I_app_m':1, #uA/cm2
    'I_app_a':1, #uA/cm2
    'cof':500,#Hz
    'C':1,#uF/cm2
    #Conductances
    'g_Na':100, #mS/cm2
    'g_K':200, #mS/cm2
    'g_L':0.1, #mS/cm2
    #Activation and inactivation parameters
    #I_Na
    'alpha_mV':54,
    'beta_mV':27,
    'alpha_hV':50,
    'beta_hV':27,
    # 'alpha_h_mag':0.128,
    # 'alpha_h_ts':18,
    #I_K
    'alpha_nV':52,
    'beta_nV':57,
    #Extracellular concentrations
    'Na_o':140,#mM
    'K_o':15,#mM
    #Cell properties
    'rho':4000,#check.. (cellular Area/Volume [1/cm]  Vitzthum,2009)
    #(http://nitrolab.engr.wisc.edu/pubs/vitzthumIB09.pdf)
    'remim':0.2,#Ratio ef extra to intra-cellular volume Zhang 2010,
    #( http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2920647/ )
    'P_K':0.96,
    'P_Na':0.04,
    'F':96484.6,#C/mol (Faraday constant)
    'R':8314.4,#mJ/[mol*K] (Ideal gas constant)
    'T':273.15+20,#K Temperature
    #Pump parameters
    'I_max':40,#uA/uF... Paper Luo y Rudi 1.5!!
    #Simple pump parameters
    'Na_sens0':20,
    'Na_sens':0.1,
    'parstNa':3,#Pump stoichiometry
    'parstK':2,#Pump stoichiometry
    #Relative Timing of ionic accumulation
    'r_fast_vs_ion':0.001,## Accounting for the fact the units of the change in concentrations is 1mM/S.. and the fast variables are in 1/ms
    #Temperature dependencies Q10s
    'T0':273.15+18,
    'Q10_gL':1.2,
    'Q10_gNa':1.2,
    'Q10_gK':1.2,
    'Q10_gPump':1.2,
    'Q10_n':2,
    'Q10_m':2,
    'Q10_h':2,
    #Initial conditions.. easy to manipulate from here at
    #the time of parameter exploration
    'm_Na0':0.1,
    'h_Na0':0.6,
    'n_K0':0.4,
    'Na_i0':10,
    'K_i0':150,#Very high...(To initialize it far from the dep block)
    'K_o0':8,
    'V0':-60
    }

# Mile traub model with simple pump[Na sensitive] definition
class neuron_MTM_W_sPNaS_sICD():

        def __init__(self, p):
            self.s_model_tag="Traub-Miles Model with Na sensitive Pump, and temperature dependencies"
            self.s_model_reference=""
            self.p = p
            self.n_state_vars=7
            self.s_state_vars=["m_Na", "h_Na", "n_K", "Na_i",
            "K_i", "K_o", "V"]
            self.s_currents=["I_Na", "I_K", "I_1", "i_p", "I_L","I_L_K", "I_L_Na"]
            self.s_concentrations=["Na_i","K_i","K_o"]
            self.s_resting_membrane_potentials=["E_Na","E_K","E_L"]
            self.step_protocol=[]
            self.current_state=[self.p['m_Na0'],self.p['h_Na0'],self.p['n_K0'],
            self.p['Na_i0'],self.p['K_i0'],self.p['K_o0'],self.p['V0']]
            self.ode_solver="odeint"
            self.noisy_current=[]
        def resting_membrane_potentials(self,v_State_Vars):
            #Resting potential for each neuron
            m_Na, h_Na, n_K, Na_i, K_i, K_o, V = v_State_Vars
            E_Na=self.p['R']*self.p['T']/self.p['F']*np.log(self.p['Na_o']/Na_i)
            E_K=self.p['R']*self.p['T']/self.p['F']*np.log(K_o/K_i)
            E_L=self.p['R']*self.p['T']/self.p['F']*np.log((K_o*self.p['P_K']+
            self.p['Na_o']*self.p['P_Na'])/(K_i*self.p['P_K']+Na_i*self.p['P_Na']))
            return E_Na, E_K, E_L

        def neuron_currents(self,v_State_Vars,t=[],i_p_f=[]):
            #Resting potential for each neuron
            m_Na, h_Na, n_K, Na_i, K_i, K_o, V = v_State_Vars# Voltage it's always the last state variable
            E_Na, E_K, E_L=self.resting_membrane_potentials(v_State_Vars)
            #Pump
            if i_p_f==[]:
                if Na_i<=self.p['Na_sens0']/2:
                    i_p=0
                if Na_i>self.p['Na_sens0']/2:
                    temp_factor_p=self.p['Q10_gPump']**(
                    (self.p['T']-self.p['T0'])/10)
                    i_p=self.p['I_max']*(1/(1+np.exp(-self.p['Na_sens']*(
                    Na_i-self.p['Na_sens0'])))-1/(1+np.exp(-self.p['Na_sens']*(
                    -self.p['Na_sens0']/2))))*temp_factor_p
            else:
                i_p=i_p_f(t) ##Argument to force the pump to be f_i_p value
            # Fast Na current
            I_Na=self.p['g_Na']*m_Na**3*h_Na*(V-E_Na)*self.p['Q10_gNa']**(
            (self.p['T']-self.p['T0'])/10)
            #K+ delayed rectifier
            I_K=self.p['g_K']*n_K**4*(V-E_K)*self.p['Q10_gK']**(
            (self.p['T']-self.p['T0'])/10)
            #Leak
            I_L_K=self.p['g_L']*self.p['P_K']*(V-E_K)*self.p['Q10_gL']**(
            (self.p['T']-self.p['T0'])/10)#[uA/cm2]
            I_L_Na=self.p['g_L']*self.p['P_Na']*(V-E_Na)*self.p['Q10_gL']**(
            (self.p['T']-self.p['T0'])/10)#[uA/cm2]
            I_L=I_L_K+I_L_Na
            I_1=I_L+i_p*(self.p['parstNa']-self.p['parstK'])#[uA/cm2]
            return I_Na, I_K, I_1, i_p, I_L, I_L_K, I_L_Na

        def activation_inactivation_Dynamics(self,v_State_Vars):
            m_Na, h_Na, n_K, Na_i, K_i, K_o, V = v_State_Vars
            #Activation and inactivation variables
            #fast Na current
            alpha_m=0.32*(V+self.p['alpha_mV'])/(1-np.exp(
            -(V+self.p['alpha_mV'])/4))*self.p['Q10_m']**(
            (self.p['T']-self.p['T0'])/10)
            beta_m=0.28*(V+self.p['beta_mV'])/(np.exp(
            (V+self.p['beta_mV'])/5)-1)*self.p['Q10_m']**(
            (self.p['T']-self.p['T0'])/10)
            alpha_h=0.128*np.exp(-(V+self.p['alpha_hV'])/18)*self.p['Q10_h']**(
            (self.p['T']-self.p['T0'])/10)
            beta_h=4/(1+np.exp(-(V+self.p['beta_hV'])/5))*self.p['Q10_h']**(
            (self.p['T']-self.p['T0'])/10)
            #K+ delayed rectifier
            alpha_n=0.032*(V+self.p['alpha_nV'])/(1
            -np.exp(-(V+self.p['alpha_nV'])/5))*self.p['Q10_n']**(
            (self.p['T']-self.p['T0'])/10)
            beta_n=0.5*np.exp(-(V+self.p['beta_nV'])/40)*self.p['Q10_n']**(
            (self.p['T']-self.p['T0'])/10)
            return alpha_m, beta_m, alpha_h, beta_h, alpha_n, beta_n

        def neuron_changing_state(self,v_State_Vars,t,I_app,i_p_f):
            m_Na, h_Na, n_K, Na_i, K_i, K_o, V = v_State_Vars
            I_Na, I_K, I_1, i_p, I_L, I_L_K, I_L_Na=self.neuron_currents(v_State_Vars,t,i_p_f)
            alpha_m,beta_m,alpha_h,beta_h,alpha_n,beta_n=self.activation_inactivation_Dynamics(v_State_Vars)
            #Change of state variables
            dm_Na=alpha_m*(1-m_Na)-beta_m*m_Na
            dh_Na=alpha_h*(1-h_Na)-beta_h*h_Na
            dn_K=alpha_n*(1-n_K)-beta_n*n_K
            dNa_i=(self.p['rho']/self.p['F']*(
            -self.p['parstNa']*i_p-I_Na-I_L_Na)
            )*self.p['r_fast_vs_ion']
            dK_i=(self.p['rho']/self.p['F']*(self.p['parstK']*i_p
            -I_K-I_L_K))*self.p['r_fast_vs_ion']#Shut off o on..
            dK_o=-(dK_i)*self.p['remim']#Shut off o on..
            dV=1/self.p['C']*(I_app(t)-I_1-I_Na-I_K)
            return dm_Na, dh_Na, dn_K, dNa_i, dK_i, dK_o, dV

        def stimulate_neuron(self,t,c_ini,I_stim,i_p_f=[]):
            import copy
            sol= odeint(self.neuron_changing_state, c_ini,
            t, args=(I_stim,i_p_f))
            s_results=[]
            s_results=copy.copy(self.s_state_vars)
            s_results.append("t")
            t_prov=t[...,None]
            v_results=np.append(sol,t_prov,1)
            return s_results, v_results