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