############################################### # simulation code for Fig. 7; # # with "chi=0.2" also Fig. 8b can be produced # ######################################################## # t: time in [msec] # # v: membrane potential in [mV] # # n,h: gating variables in [1] # # n_ki: amount of potassium in the ICS in [fmol] # # n_cli: amount of chloride in the ICS in [fmol] # # dnk: amount of potassium that is # # exchanged with glia (negative # # values when potassium goes into # # glia cell) in [fmol] # # vli,vlg: ICS and glia volume in [um^3] # ######################################################## ################################################## # rate equations used by solver # # factor 1000. converts from [x/msec] to [x/sec] # ################################################## v' = 1000. * V_DOT n' = 1000. * N_DOT h' = 1000. * H_DOT n_ki' = 1000. * N_KI_DOT n_cli' = 1000. * N_CLI_DOT dnk' = 1000. * DNK_DOT vli' = 1000. * VLI_DOT vlg' = 1000. * VLG_DOT ############################################ # initial conditions: normal resting state # ############################################ init v=-67. init n=0.070 init h=0.978 init n_ki=277.7 init n_cli=21.7 init dnk=0 init vli=2160. init vlg=2160. ################################################################### # by the parameter delta we set the duration of pump interruption # # and interruption of glial ion regulation in [sec]; # # f_pmp/gl: pump/glia factor in [1] # ################################################################### par delta=20 f_pmp = (heav(50-t) + heav(t-50-delta)) f_gl = (heav(50-t) + heav(t-50-delta)) ################################################################################ # chi sets the fraction of chloride cotransport with glial potassium buffering # ################################################################################ par chi=0.8 ######################################## # Hodgkin-Huxley like gating functions # # adiabatic value for "M" # ######################################## AN = 0.01 * (v + 34.0) / (1.0 - exp(-0.1 * (v + 34.0))) BN = 0.125 * exp(-(v + 44.0) / 80.0) AM = 0.1 * (v + 30.0) / (1.0 - exp(-0.1 * (v + 30.0))) BM = 4.0 * exp(-(v + 55.0) / 18.0) AH = 0.07 * exp(-(v + 44.0) / 20) BH = 1.0 / (1.0 + exp(-0.1 * (v + 14.0))) M = AM / (AM + BM) ##################################################################### # - resting state ion and impermeant particle amounts in [fmol] # # in the ICS (...i0) and the ECS (...e0) # # - "k,na,cl,imp" means potassium, sodium, chloride and impermeant # ##################################################################### n_ki0 = 277.7 n_ke0 = 2.8 n_nai0 = 54.6 n_nae0 = 91.3 n_cli0 = 21.7 n_cle0 = 89.8 n_impi = 318. n_impe = 40. ############################################################################ # electroneutrality (first line) and mass conservation (second to fourth) # # conditions to compute ion amounts other than intracellular potassium and # # chloride # ############################################################################ n_nai = n_nai0 + n_ki0 - n_ki - n_cli0 + n_cli dnna =-dnk * (1-chi) dncl = dnk * chi n_nae = n_nae0 + n_nai0 - n_nai + dnna n_ke = n_ke0 + n_ki0 - n_ki + dnk n_cle = n_cle0 + n_cli0 - n_cli + dncl ################################################################################## # ECS volume "vle" in [um^3]: # # "vle_osm" is the value which follows from the other volumes and from # # conservation of the initial total volume "vl_tot0" - a value based only on # # osmosis. "vle" equals this value as long as "vle_osm" is not too small. For # # very small values of "vle_osm" the below function ensures that "vle > vle_osm" # # note: the equation for "vle" can be converted to Eq. (49) in the manuscript by # # substituting "vle_osm" with "vl_tot0*n_e/n_tot" # ################################################################################## # n_i/e/g/tot: amount of particles in ICS/ECS/glia/whole system # # vl_tot0: initial total volume # ################################################################# n_i = n_nai + n_ki + n_cli + n_impi n_e = n_nae + n_ke + n_cle + n_impe n_g = (n_ki0 + n_nai0 + n_cli0 + n_impi) - (dnk + dnna + dncl) n_tot = n_i + n_e + n_g vl_tot0 = 5040. vle_osm = vl_tot0 - vli - vlg p1 = 0.2 p2 = 5 p3 = 0.93 p4 = 0.2 p5 = 0.21 p6 =-0.095 vle = 1000. * ((p3 * (vle_osm/1000.-p6) - p4) / (1 + exp((p1-(vle_osm/1000.-p6)) * p2)) + p5) vli_inf = n_i * vle / n_e vlg_inf = n_g * vle / n_e ################################################# # equilibrium volume "vli/g_inf" for glia and # # neuron based on osmotic balance of these # # compartments with the ECS # ################################################# ######################################## # ion concentrations in [mM]=[mMol/l] # ######################################## nai = n_nai / vli * 1000. nae = n_nae / vle * 1000. ki = n_ki / vli * 1000. ke = n_ke / vle * 1000. cli = n_cli / vli * 1000. cle = n_cle / vle * 1000. ##################### # Nernst potentials # ##################### EK = 26.64 * log(ke / ki) ENA = 26.64 * log(nae / nai) ECL =-26.64 * log(cle / cli) ############################################################################## # - different leak and gated currents "IION_l/g" in [umA/cm^2] # # - leak and gated conductances for each channel "gion_l/g" in [mS/cm^2] # # - Na/K-exchange pump current "IP" with maximal turnover rate "max_p" # ############################################################################## gna_l = 0.0175 gna_g = 100. gk_l = 0.05 gk_g = 40. gcl_l = 0.05 max_p = 6.8 INA_l = gna_l * (v - ENA) INA_g = gna_g * M**3 * h * (v - ENA) IK_l = gk_l * (v - EK) IK_g = gk_g * n**4 * (v - EK) ICL_l = gcl_l * (v - ECL) IP = max_p / (1.0 + exp((25 - nai)/3.)) / (1. + exp(5.5 - ke)) * f_pmp ##################################### # full sodium and potassium current # ##################################### INA = INA_l + INA_g + 3. * IP IK = IK_l + IK_g - 2. * IP ################### # glial buffering # ################### k1 = 1.75e-3 k2 = 6.2e-4 J_gl = (k2 - k1 / (1.0 + exp((5.5-ke)/2.5))) * f_gl ############################# # list of all changes rates # ##################################################################### # c: membrane capacitance in [uF/cm^2] # # conv: conversion from currents to fluxes in [fmol/msec*cm^2/uA] # # phi: gating timescale parameter in [1/msec] # # t_vl: timescale of volume dynamics in [msec] # ##################################################################### c = 1 conv = 9.55589e-05 phi = 3 t_vl = 250 V_DOT = -1. / c * (INA + IK + ICL_l) N_DOT = phi * (AN * (1 - n) - BN * n) H_DOT = phi * (AH * (1 - h) - BH * h) N_KI_DOT = -CONV * IK N_CLI_DOT = CONV * ICL_l DNK_DOT = J_gl VLI_DOT = 1. / t_vl * (vli_inf - vli) VLG_DOT = 1. / t_vl * (vlg_inf - vlg) #################################### # auxiliary variables for plotting # #################################### aux _ki = ki aux _ke = ke aux _nai = nai aux _nae = nae aux _cli = cli aux _cle = cle aux _EK = EK aux _ENA = ENA aux _ECL = ECL vli0 = 2160. vlg0 = 2160. vle0 = vl_tot0 - vli0 - vlg0 aux vli_line = vli aux vle_line = (vli+vle) aux vlg_line = (vli+vle+vlg) aux dvli = (vli - vli0) / vli0 * 100 aux dvle = (vle - vle0) / vle0 * 100 aux dvlg = (vlg - vlg0) / vlg0 * 100 ############ # numerics # ############ @ meth=cvode @ dt=1e-4 @ maxstor=10000000, bounds=10000000 @ total=500 @ bell=0 ############################################# # plot options corresponding to Fig. 7a/b/c # ############################################# @ xhi=500 @ nplot=4, yp1=v, yp2=_EK, yp3=_ENA, yp4=_ECL, ylo=-150, yhi=160 #@ nplot=3, yp1=vli_line, yp2=vle_line, yp3=vlg_line, ylo=0, yhi=6000 #@ nplot=3, yp1=dvli, yp2=dvle, yp3=dvlg, ylo=-90, yhi=50 done