############################################################################################### # Code for fixed point continuations/bifurcation diagrams used in # # "Large Extracellular Space Leads to [...] Ischemic Injury [...]" by Hubel, Ullah and Andrew # ############################################################################################### ############################################################################################################## ############################################################################################################## ### ### ### The fixed point curves from Figs. 4 and 5 can be obtained as follows: ### ### ### ### 1.) open file with XPPAUT ### ### 2.) run simulation twice to make sure the system is in its fixed point: ### ### click "Initialconds" + "(G)o"; "Initialconds" + "(L)ast" ### ### (keyboard shortcuts: "I" + "G", "I" + "L") ### ### 3.) open AUTO interface: ### ### click "File" + "Auto" ### ### (keyboard shortcut: "F" + "A") ### ### 4.) run 'forward' fixed point continuation: ### ### click "Run" + "Steady state" ### ### (keyboard shortcut: "R" + "S") ### ### 5.) grab point to start 'backward' continuation if desired: ### ### click "Grab", then navigate along the curve with "tab" key, press "enter" to choose point ### ### 6.) set continuation step size to negative: ### ### click "Numerics" and change 'Ds:0.002' to 'Ds:-0.002', click "Ok" ### ### 7.) run backward coninuation by clicking "Run" ### ### 8.) save the fixed point curves and bifurcation information: ### ### click "File" + "All info", choose a filename and click "Ok" ### ### ### ### Remark for Fig. 4: ### ### Start from the 'physiological' (i.e. upper set of) initial conditions and 'par vle=720'. ### ### Change 'PARMIN=-250' to 'PARMIN=-50' and accordingly 'XAUTOMIN=-250' to 'XAUTOMIN=-50'. ### ### For more negative parameter values the continuation produces convergence errors "MX" at ### ### extreme negative values of 'V' ### ### ### ### Remark for Fig. 5a: ### ### Start from the 'FES' (i.e. lower set of) initial conditions and 'par vle=3700' to obntain ### ### the upper branch of the continuation. Use 'physiological' initial conditions and 'par vle=3700' ### ### to obtain the lower loop of the bifurcation diagram ### ### ### ############################################################################################################## ############################################################################################################## #################################################################################### # membrane potential 'V' in [mV] # # gating variables 'n/h' in [1] # # intracellular ion concentrations 'ki,cli' in [mM=mMol/l] # # rates of change 'X_DOT' in [.../msec] # # (factor 1000. converts to seconds) # #################################################################################### V' = 1000. * V_DOT n' = 1000. * N_DOT h' = 1000. * H_DOT ki' = 1000. * KI_DOT cli' = 1000. * CLI_DOT ###################### # Initial conditions # ###################### # physiological # ################# # init v=-67.056664 # init n=0.070174225 # init h=0.97820824 # init ki=128.56937 # init cli=10.061391 ###################### # FES for 'vle=3700' # ###################### init v=-20.288656 init n=0.66411209 init h=0.057942949 init ki=77.408791 init cli=47.911449 ######################################################################### # MAIN BIFURCATION PARAMETER: # # amount of potassium exchanged with external reservoir 'dnk' in [fmol] # ######################################################################### par dnk=0 ######################################## # Extracellular volume 'vle' in [um^3] # ######################################## # par vle=720 par vle=3700 ################################## # Pump strength 'max_p' constant # ################################## max_p = 6.8 ################################### # Hodgkin-Huxley gating functions # ################################### 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) ##################################### # ion concentrations in [mM=mMol/l] # ##################################################################### # intracellular sodium 'NAI' from electroneutrality # # extracellular concentrations 'NAE,KE,CLE' from mass conservation # # normal resting values 'ki0,ke0 ...' given # ##################################################################### vli = 2160 ki0 = 128.56935 ke0 = 3.9962559 nai0 = 25.279156 nae0 = 126.84917 cli0 = 10.055541 cle0 = 124.71021 NAI = nai0 + ki0 - ki - cli0 + cli NAE = nae0 + (nai0 - nai) * vli/vle KE = ke0 + (ki0 - ki ) * vli/vle + dnk/vle * 1e3 CLE = cle0 + (cli0 - cli) * vli/vle ############################# # Nernst potentials in [mV] # ############################# EK = 26.64 * log(ke / ki) ENA = 26.64 * log(nae / nai) ECL =-26.64 * log(cle / cli) ############################################################################ # different types of 'l'eak and 'g'ated currents 'I(ION)_l/g' in [uA/cm^2] # # different channel conductances 'g(ion)_l/g' in [mS/cm^2] # # Na/K-exchange pump current 'IP' in [uA/cm^2] # ############################################################################ gna_l = 0.0175 gna_g = 100. gk_l = 0.05 gk_g = 40. gcl_l = 0.05 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)) INA = INA_l + INA_g + 3. * IP IK = IK_l + IK_g - 2. * IP ############################# # Full list of change rates # ########################################################## # membrane capacitance 'C' in [uF/cm^2] # # conversion factor 'conv' in [XXX] # # conventional time scale parameter 'phi' in [1/msec] # ########################################################## c = 1 conv = 9.55589e-2 phi = 3 V_DOT = -1. / c * (INA + IK + ICL_l) N_DOT = phi * (AN * (1 - n) - BN * n) H_DOT = phi * (AH * (1 - h) - BH * h) KI_DOT = -CONV/vli * IK CLI_DOT = CONV/vli * ICL_l #################################################### # auxiliary variables for data output and 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 #################################### # Numerical parameters: simulation # #################################### @ meth=stiff @ dt=5e-3 @ maxstor=10000000, bounds=10000000 @ total=200 @ bell=0 ############################################################# # Parameters for fixed point continuation in AUTO interface # ############################################################# @ NTST=50, NMAX=115000, NPR=115000 @ DS=0.002, DSMIN=0.001, DSMAX=0.005 @ PARMIN=-250, PARMAX=50 @ AUTOXMIN=-250, AUTOXMAX=50, AUTOYMIN=-150, AUTOYMAX=50 @ EPSL=0.001, EPSU=0.001, EPSS=0.001 ################# # Plot settings # ########################################################### # '_MAX_P' as a guide to the eye to see pump interruption # ########################################################### @ xhi=200 @ nplot=3, yp1=v, yp2=_EK, yp3=_ENA, ylo=-150, yhi=160 done