######################################################################################## # This code is archived with “The Role of Glutamate in Neuronal Ion Homeostasis: A Case Study of Spreading # Depolarization” by Hubel et al. PLoS Computational Biology (2017). # The software used is AUTO. # Feel free to use/modify this code for your own research. Please, cite our paper if you use this code. ######################################################################################## # RATE EQUATIONS FOR SD WITH GLUTAMATE: # # 'v' membrane potential in [mV] # # 'n/h_HH' Hodgkin-Huxley gating variables in [1] # # 'r_NMDA/AMPA' NMDA/AMPA gating variables in [1] # # 'nki/ncli' amount of K/Cl in ICS in [fmol] # # 'dnk_b/g' amount of K that is exchanged with external bath/glia in [fmol] # # 'ngl_c/e' amount of glutamate in all clefts/ECS in [fmol] # ######################################################################################## v' = 1000. * v_DOT n_HH' = 1000. * n_HH_DOT h_HH' = 1000. * h_HH_DOT r_AMPA' = 1000. * r_AMPA_DOT * f_receptor r_NMDA' = 1000. * r_NMDA_DOT * f_receptor nki' = 1000. * nki_DOT ncli' = 1000. * ncli_DOT dnk_b' = 1000. * dnk_b_DOT dnk_g' = 1000. * dnk_g_DOT ngl_c' = 1000. * ngl_c_DOT ngl_e' = 1000. * ngl_e_DOT init v=-71.1 init n_HH=0.07143807 init h_HH=0.9774849 init r_AMPA=0. init r_NMDA=0. init nki=1050. init ncli=67.5 init dnk_b=0. init dnk_g=0. init ngl_c=0. init ngl_e=0. k_rep = 0.001 ngl_rel' = 1000. * (J_rel - k_rep / ngl_tot * (ngl_c_i + ngl_e_i + ngl_c_g + ngl_e_g)) ngl_c_i' = 1000. * (V_c_i - k_rep / ngl_tot * ngl_c_i) ngl_c_g' = 1000. * (V_c_g - k_rep / ngl_tot * ngl_c_g) ngl_c_e' = 1000. * Jglu_diff ngl_e_i' = 1000. * (V_e_i - k_rep / ngl_tot * ngl_e_i) ngl_e_g' = 1000. * (V_e_g - k_rep / ngl_tot * ngl_e_g) init ngl_rel=0 init ngl_c_i=0 init ngl_c_g=0 init ngl_c_e=0 init ngl_e_i=0 init ngl_e_g=0 ############################################################################# # GLUTAMATE UPTAKE/RELEASE PARAMETERS: # # 'exp_rel' exponent for 'v'-dependence of glutamate release in [1] # # 'Vmax' maximal velocity of glutamate uptake by neuron in [mM/msec] # ############################################################################# par f_upt=1 par f_co=1 par Vmax=0.03 exp_rel = 2 ################################ # switch for receptor dynamics # ################################ f_receptor = 1 ####################################### # SD from extrac. potassium elevation # ####################################### f_OGD = 1 par KB=15. ############### # SD from OGD # ############### # par KB=4. # par delta_t=15 # par t_OGD=40 # # t1 = 20 # t2 = t1 + delta_t # t3 = t2 + t_OGD # t4 = t3 + delta_t # f01_OGD = 1 # f12_OGD = heav(t - t1) * (-(t - t1) / (t2 - t1)) # f23_OGD = heav(t - t2) * (t - t2) / (t2 - t1) # f34_OGD = heav(t - t3) * (t - t3) / (t4 - t3) # f4x_OGD = heav(t - t4) * (-(t - t4) / (t4 - t3)) # f_OGD = f01_OGD + f12_OGD + f23_OGD + f34_OGD + f4x_OGD ############################################### # GATING FUNCTIONS: # # Hodgkin-Huxley (HH) and NMDA/AMPA receptors # # (gating time constants in [msec]) # ############################################### an_HH = 0.01 * (v + 34.0) / (1.0 - exp(-0.1 * (v + 34.0))) bn_HH = 0.125 * exp(-(v + 44.0)/80.0) am_HH = 0.1 * (v + 30.0) / (1.0 - exp(-0.1 * (v + 30.0))) bm_HH = 4.0 * exp(-(v + 55.0)/18.0) ah_HH = 0.07 * exp(-(v + 44.0)/20.0) bh_HH = 1.0 / (1.0 + exp(-0.1 * (v + 14.0))) m_HH = am_HH / (am_HH + bm_HH) ar_AMPA = 1.1 br_AMPA = 0.19 ar_NMDA = 0.072 br_NMDA = 0.0066 ############################################################################################# # MORPHOLOGY: # # 'vl_i0/e0/g0' initial/resting volume of ICS, ECS and glia in [um^3] # # 'A_m' membrane surface area in [um^2] (same for neuron and glia) # # 'h','r' height and radius of synaptic cleft in [um] # # 'vl_en' volume inside the glial envelope in [um^3] # # 'A_sig' cross section area for glutamate diffusion from cleft into ECS in [um^2] # ############################################################################################# vl_i0 = 7500. vl_e0 = 2500. vl_g0 = 7500. vl_tot = vl_i0 + vl_e0 + vl_g0 A_m = 18000. r = 100e-3 h = 20e-3 vl_en = 6 * pi * r**2 * h A_sig = 4 * pi * r**2 * 0.05 ################################################################################## # ION CONENTRATIONS: # # 'ioni/e0' initial/resting conc. in ICS/ECS in [mM=mmol/l] # # 'xi' conc. of impermeant uncharged matter in ICS in [mM] # # 'nioni/e0' initial/resting ion amount in ICS/ECS in [fmol] # # 'nani/e' amount of impermeant anions in ICS/ECS in [fmol] # # 'nxi/e' amount of impermeant uncharged particles in ICS/ECS in [fmol] # # 'ni0/e0/g0' initial amount of all particles in ICS/ECS/glia in [fmol] # # 'dnna_b/g' amount of Na that is exchanged with external bath/glia in [fmol] # # 'dncl_b/g' amount of Cl that is exchanged with external bath/glia in [fmol] # # 'NNAI/E' current amount of Na in ICS/ECS in [fmol] # # 'NKE/CLE' current amount of K and Cl in ECS in [fmol] # # 'ni/e/g' current amount of all particles in ICS/ECS/glia in [fmol] # # 'n_tot' total amount of particles in the whole system in [fmol] # # 'vl_i/e/g' current volume of ICS/ECS/glia in [um^3] # # 'IONI/E' current ion conc. in ICS/ECS in [mM=mmol/l] # # (conversion factor 1e3 from [fmol/um^3] to [mM] for 'KE') # ################################################################################## nai0 = 15. nae0 = 144. ki0 = 140. ke0 = 4. cli0 = 9. cle0 = 130. xi = 5. nki0 = ki0 * vl_i0 * 1e-3 nke0 = ke0 * vl_e0 * 1e-3 nnai0 = nai0 * vl_i0 * 1e-3 nnae0 = nae0 * vl_e0 * 1e-3 ncli0 = cli0 * vl_i0 * 1e-3 ncle0 = cle0 * vl_e0 * 1e-3 nani = nki0 + nnai0 - ncli0 nane = nke0 + nnae0 - ncle0 nxi = xi * vl_i0 * 1e-3 nxe = vl_e0/vl_i0 * (nxi + nani + nki0 + nnai0 + ncli0) - (nane + nke0 + nnae0 + ncle0) ni0 = nnai0 + nki0 + ncli0 + nani + nxi ne0 = nnae0 + nke0 + ncle0 + nane + nxe ng0 = ni0 / vl_i0 * vl_g0 dncl_g = 0.8 * dnk_g dnna_g =-0.2 * dnk_g NNAI = nnai0 + nki0 - nki - ncli0 + ncli NNAE = nnae0 + nnai0 - NNAI - dnna_g NKE = nke0 + nki0 - nki - dnk_g - dnk_b NCLE = ncle0 + ncli0 - ncli - dncl_g ni = NNAI + nki + ncli + nani + nxi ne = NNAE + NKE + NCLE + nane + nxe ng = ng0 + dnna_g + dnk_g + dncl_g n_tot = ni + ne + ng vl_i = ni / n_tot * vl_tot vl_e = ne / n_tot * vl_tot vl_g = ng / n_tot * vl_tot NAI = NNAI / vl_i * 1e3 KI = NKI / vl_i * 1e3 CLI = NCLI / vl_i * 1e3 NAE = NNAE / vl_e * 1e3 KE = NKE / vl_e * 1e3 CLE = NCLE / vl_e * 1e3 ####################################################### # NERNST POTENTIALS: # # 'EK/NA/CL' Nernst potetnials for K/Na/Cl in [mV] # ####################################################### EK = 26.64 * log(KE / KI) ENA = 26.64 * log(NAE / NAI) ECL =-26.64 * log(CLE / CLI) ####################################################################################### # CURRENTS ON NEURAL MEMBRANE: # # 'gk/na_g' maximal conductances of gated HH-like K/Na channels in [mS/cm^2] # # 'gk/na/cl_l' conductances of K/Na/Cl leak channels in [mS/cm^2] # # 'nsyn_AP/SD' number of synapses activated in action potential (AP)/SD in [1] # # 'g_NMDA/AMPA_AP' effective NMDA/AMPA receptor conductance in single AP in [mS] # # 'g_NMDA/AMPA' effective NMDA/AMPA conductance density for SD in [mS/cm^2] # # 'max_p' maximal pump rate for K/Na-exchange in [uA/cm^2] # # 'mg' Mg concentration in ECS in [mM] # # 'f_NMDA' factor for voltage-dependent Mg block in NMDA channels in [1] # # 'IK/NA_g' gated K/Na current in [uA/cm^2] # # 'IK/NA/CL_l' leak K/Na/Cl current in [uA/cm^2] # # 'IK/NA_NMDA' K/Na current through NMDA receptor gates in [uA/cm^2] # # 'IK/NA_AMPA' K/Na current through AMPA receptor gates in [uA/cm^2] # # 'IP' K/Na-exchange pump current in [uA/cm^2] # # 'INA/K' sum of all Na/K currents in [uA/cm^2] # ####################################################################################### gk_g = 40. gna_g = 100. gna_l = 0.0135 gk_l = 0.05 gcl_l = 0.05 nsyn_AP = 20 nsyn_SD = 5000 par g_NMDA_AP=1e-7 par g_AMPA_AP=3.5e-7 g_AMPA = g_AMPA_AP / (A_m * 1e-8) * nsyn_SD / nsyn_AP g_NMDA = g_NMDA_AP / (A_m * 1e-8) * nsyn_SD / nsyn_AP mg = 1.2 f_NMDA = 1.0 / (1.0 + 0.33 * mg * exp(-(0.07*v + 0.7))) max_p = 6.46 INA_g = gna_g * m_HH**3 * h_HH * (v-ENA) IK_g = gk_g * n_HH**4 * (v-EK) INA_l = gna_l * (v-ENA) IK_l = gk_l * (v-EK) ICL_l = gcl_l * (v-ECL) IK_NMDA = g_NMDA * r_NMDA * (v-EK) * f_NMDA INA_NMDA = g_NMDA * r_NMDA * (v-ENA) * f_NMDA IK_AMPA = g_AMPA * r_AMPA * (v-EK) INA_AMPA = g_AMPA * r_AMPA * (v-ENA) IP = max_p / (1.0 + exp((15 - nai)/3.)) / (1. + exp(5.5 - ke)) * f_OGD INA = INA_l + INA_g + INA_AMPA + INA_NMDA + 3. * IP IK = IK_l + IK_g + IK_AMPA + IK_NMDA - 2. * IP ############################################################################################# # GLUTAMATE RELEASE IN ONE SYNAPSE # # 'ngl_tot/i' total/intracellular amount of glutamate in [fmol] # # 'v_cr' critical mebrane potential for glutamate release in [mV] # # 'v_hi' highest value of membrane potential seen in the model in [mV] # # 'rel_max' maximal (i.e. for 'v=v_hi') release rate in [fmol/msec] # # 'f_rel' factor for 'V'-dependence of glutamate release in [1] # # 'J_rel' glutamate release flux into all synaptic clefts involved in SD [fmol/msec] # ############################################################################################# rel_max = 1.5e-5 par ngl_tot=10 ngl_avl = ngl_tot - ngl_rel v_cr =-50 v_hi = 50 f_rel = heav(v-v_cr) * ((v-v_cr)/(v_hi-v_cr))**exp_rel J_rel = nsyn_SD * rel_max * f_rel * ngl_avl / ngl_tot #################################################################################### # GLUTAMATE UPTAKE: # # 'gl_c/e' glutamate concentration in clefts/ECS in [mM] # # 'k_m' affinity pf uptake system in [mM] # # 'Vmax_c/e_i' maximal uptake velocity from cleft/ECS into ICS in [mM/msec] # # 'Vmax_c/e_g' maximal uptake velocity from cleft/ECS into glia in [mM/msec] # # 'F' Faraday's constant in [A*sec/mol] # # 'conv' current-to-flux conversion s.th. 'conv/vl_i*INA' is in [mM/msec] # # 'nsyn_AP/SD' number of synapses activated in action potential (AP)/SD in [1] # # 'V_c/e_i' uptake velocity from cleft/ECS into ICS in [fmol/msec] # # 'V_c/e_g' uptake velocity from cleft/ECS into glia in [fmol/msec] # # 'JNa/Cl/K_co' glutamate cotransport flux of Na/Cl/K into ICS in [mM/msec] # # 'INa/Cl/K_co' cotransport Na/Cl/K currents on neural membrane in [uA/cm^2] # #################################################################################### gl_c = ngl_c / (vl_en * nsyn_SD) * 1e3 gl_e = ngl_e / vl_e * 1e3 k_m = 3e-2 Vmax_c_i = f_upt * Vmax Vmax_e_i = f_upt * Vmax * 0.12 * vl_e0 / vl_e Vmax_c_g = f_upt * Vmax * 4 Vmax_e_g = f_upt * Vmax * 4 * 0.24 * vl_e0 / vl_e F = 96485 conv = A_m / F * 10 V_c_i = Vmax_c_i * gl_c / (gl_c + k_m) * vl_en * nsyn_SD * 1e-3 * f_OGD V_e_i = Vmax_e_i * gl_e / (gl_e + k_m) * vl_e * 1e-3 * f_OGD V_c_g = Vmax_c_g * gl_c / (gl_c + k_m) * vl_en * nsyn_SD * 1e-3 * f_OGD V_e_g = Vmax_e_g * gl_e / (gl_e + k_m) * vl_e * 1e-3 * f_OGD JNa_co = (V_c_i + V_e_i) * 3 * f_co JCl_co = (V_c_i + V_e_i) * f_co JK_co =-(V_c_i + V_e_i) * f_co INa_co = JNa_co / conv * 1e3 ICl_co =-JCl_co / conv * 1e3 IK_co = JK_co / conv * 1e3 ################################################################################ # DIFFUSION PROCESSES: # # 'lambda' strength of bath coupling in [1/msec] # # 'kbath' level of K concentration of bath in [mM] # # 'JK_diff' diffusive K flux between ECS and bath in [mM/msec] # # 'D_glu' glutamate diffusion constatn in [um^2/msec] # # 'dx' distance between cleft and staionary bath concentraion in [um] # # 'Jglu_diff' diffusive glutamate flux from cleft to ECS in [fmol/msec] # ################################################################################ par lambda=0.0001 JK_diff =-lambda * (KB - KE) * f_OGD D_glu = 0.3 dx = 20 Jglu_diff =-D_glu * nsyn_SD / dx * A_sig * (ngl_e/vl_e - ngl_c/(vl_en*nsyn_SD)) ################################################################### # GLIAL BUFFERING: # # 'k1' buffering constant in [mM/msec] # # 'k2' backward buffering rate in [mM/msec] # # 'dnk_max' maximal amount of K that can go into glia in [fmol] # # 'JK_glia' K uptake flux in [fmol/msec] # ################################################################### k1 = 1.44e-2 k2 = 5.1e-3 dnk_max = 350 JK_glia =-(k2 - k1 / (1.0 + exp((5.5-KE)/2.5)) * (dnk_max - dnk_g)/dnk_max) * vl_e0 * 1e-3 * f_OGD ########################################################## # RATE FUNCTIONS USED BY SOLVER: # # 'C_m' capacitance of neural membrane in [uF/cm^2] # # 'phi' conventional time scale parameter in [1/msec] # ########################################################## C_m = 1 phi = 3 v_DOT =-1 / C_m * (INA + IK + ICL_l + INA_co + IK_co + ICL_co) n_HH_DOT = phi * (an_HH * (1-n_HH) - bn_HH * n_HH) h_HH_DOT = phi * (ah_HH * (1-h_HH) - bh_HH * h_HH) r_AMPA_dot = gl_c * ar_AMPA * (1.0 - r_AMPA) - br_AMPA * r_AMPA r_NMDA_dot = gl_c * ar_NMDA * (1.0 - r_NMDA) - br_NMDA * r_NMDA nki_DOT =-conv * 1e-3 * (IK + IK_co) ncli_DOT = conv * 1e-3 * (ICL_l + ICl_co) dnk_b_DOT = JK_diff * vl_e * 1e-3 dnk_g_DOT = JK_glia ngl_c_DOT = J_rel - Jglu_diff - V_c_i - V_c_g ngl_e_DOT = Jglu_diff - V_e_i - V_e_g #################################### # AUXILIARY VARIABLES FOR PLOTTING # #################################### aux _vl_i = vl_i aux _vl_e = vl_e aux _vl_g = vl_g aux _gl_c = gl_c aux _gl_e = gl_e aux _EK = EK aux _ENA = ENA aux _ECL = ECL aux _KI = ki aux _KE = KE aux _NAI = NAI aux _NAE = NAE aux _CLI = cli aux _CLE = CLE ############################## # NUMERICS AND PLOT SETTINGS # ############################## @ maxstor=10000000, bounds=10000000 @ bell=0 # @ meth=Runge-Kutta # @ dt=5e-5 #@ meth=stiff #@ dt=1e-4 @ meth=stiff @ dt=1e-3 #@ total=800 #@ xhi=800 @ total=400 @ xhi=400 #@ ylo=-300, yhi=300 #@ nplot=6 #@ yp1=v #@ yp2=_ECL #@ yp3=_CLI #@ yp4=_CLE #@ yp5=dnk_b #@ yp6=dnk_g @ ylo=-150, yhi=160 @ nplot=9 @ yp1=v @ yp2=_EK @ yp3=_ENA @ yp4=_ECL @ yp5=_vl_i @ yp6=_vl_e @ yp7=_vl_g @ yp8=_gl_c @ yp9=_gl_e #@ ylo=-0.5, yhi=200 #@ nplot=8 #@ yp1=_gl_e #@ yp2=_gl_c #@ yp3=_KI #@ yp4=_KE #@ yp5=_NAI #@ yp6=_NAE #@ yp7=_CLI #@ yp8=_CLE #@ ylo=0, yhi=10000 #@ nplot=3 #@ yp1=_vl_i #@ yp2=_vl_e #@ yp3=_vl_g done