# Size of variable arrays: sizeAlgebraic = 16 sizeStates = 5 sizeConstants = 16 from math import * from numpy import * def createLegends(): legend_states = [""] * sizeStates legend_rates = [""] * sizeStates legend_algebraic = [""] * sizeAlgebraic legend_voi = "" legend_constants = [""] * sizeConstants legend_constants[0] = "T in component parameters (kelvin)" legend_constants[1] = "K in component parameters_propagation (per_ms)" legend_constants[2] = "V_initial in component initial_conditions (mV)" legend_voi = "time in component time (ms)" legend_algebraic[6] = "INa in component INa (uA_per_cmsq)" legend_algebraic[11] = "IK in component IK (uA_per_cmsq)" legend_algebraic[13] = "Ileak in component Ileak (uA_per_cmsq)" legend_algebraic[14] = "Ii in component propagated_AP (uA_per_cmsq)" legend_algebraic[0] = "Im in component propagated_AP (uA_per_cmsq)" legend_algebraic[15] = "I in component propagated_AP (uA_per_cmsq)" legend_algebraic[12] = "g in component propagated_AP (mS_per_cmsq)" legend_algebraic[10] = "gK in component IK (mS_per_cmsq)" legend_algebraic[1] = "gNa in component INa (mS_per_cmsq)" legend_states[0] = "V in component propagated_AP (mV)" legend_algebraic[2] = "minus_V in component propagated_AP (mV)" legend_states[2] = "X in component hh_gating_variable (dimensionless)" legend_states[3] = "X in component hh_gating_variable (dimensionless)" legend_states[4] = "X in component hh_gating_variable (dimensionless)" legend_constants[3] = "n_initial in component initial_conditions (dimensionless)" legend_constants[4] = "m_initial in component initial_conditions (dimensionless)" legend_constants[5] = "h_initial in component initial_conditions (dimensionless)" legend_constants[6] = "VNa in component parameters (mV)" legend_constants[7] = "VK in component parameters (mV)" legend_constants[8] = "Vleak in component parameters (mV)" legend_constants[9] = "gNa_max in component parameters (mS_per_cmsq)" legend_constants[10] = "gK_max in component parameters (mS_per_cmsq)" legend_constants[11] = "gleak_max in component parameters (mS_per_cmsq)" legend_constants[12] = "Cm in component parameters (uF_per_cmsq)" legend_constants[13] = "V_initial in component initial_conditions (mV)" legend_algebraic[3] = "alpha in component alpha_m (per_ms)" legend_constants[14] = "phi in component temperature_factor (dimensionless)" legend_algebraic[7] = "beta in component beta_m (per_ms)" legend_algebraic[4] = "alpha in component alpha_h (per_ms)" legend_algebraic[8] = "beta in component beta_h (per_ms)" legend_algebraic[5] = "alpha in component alpha_n (per_ms)" legend_constants[15] = "phi in component temperature_factor (dimensionless)" legend_algebraic[9] = "beta in component beta_n (per_ms)" legend_states[1] = "d/dt V in component propagated_AP (mV)" legend_rates[1] = "d^2/dtV in component propagated_AP (mV)" legend_rates[3] = "d/dt X in component hh_gating_variable (dimensionless)" legend_rates[2] = "d/dt X in component hh_gating_variable (dimensionless)" legend_rates[4] = "d/dt X in component hh_gating_variable (dimensionless)" return (legend_states, legend_algebraic, legend_voi, legend_constants) def initConsts(): constants = [0.0] * sizeConstants; states = [0.0] * sizeStates; constants[0] = 18.5 constants[1] = 10.47 constants[2] = 0 states[1] = 0.0 constants[3] = 0.315 constants[4] = 0.042 constants[5] = 0.608 constants[6] = -115.0 constants[7] = 12.0 constants[8] = -10.613 constants[9] = 120.0 constants[10] = 36.0 constants[11] = 0.3 constants[12] = 1.0 constants[13] = 0.0 constants[14] = power(3.00000, (constants[0]-6.30000)/10.0000) constants[15] = power(3.00000, (constants[0]-6.30000)/10.0000) states[0] = constants[2] states[2] = constants[5] states[3] = constants[4] states[4] = constants[3] return (states, constants) def computeRates(voi, states, constants): rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic algebraic[3] = ((0.100000*(states[0]+25.0000))/(exp((states[0]+25.0000)/10.0000)-1.00000))*constants[14] algebraic[7] = 4.00000*exp(states[0]/18.0000)*constants[14] rates[3] = algebraic[3]*(1.00000-states[3])-algebraic[7]*states[3] algebraic[4] = 0.0700000*exp(states[0]/20.0000)*constants[14] algebraic[8] = (1.00000/(exp((states[0]+30.0000)/10.0000)+1.00000))*constants[14] rates[2] = algebraic[4]*(1.00000-states[2])-algebraic[8]*states[2] algebraic[5] = ((0.0100000*(states[0]+10.0000))/(exp((states[0]+10.0000)/10.0000)-1.00000))*constants[15] algebraic[9] = 0.125000*exp(states[0]/80.0000)*constants[15] rates[4] = algebraic[5]*(1.00000-states[4])-algebraic[9]*states[4] algebraic[1] = constants[9]*(power(states[3], 3.00000))*states[2] algebraic[6] = algebraic[1]*(states[0]-constants[6]) algebraic[10] = constants[10]*(power(states[4], 4.00000)) algebraic[11] = algebraic[10]*(states[0]-constants[7]) algebraic[13] = constants[11]*(states[0]-constants[8]) algebraic[14] = algebraic[6]+algebraic[11]+algebraic[13] rates[1] = constants[1]*(states[1]+(1.00000/constants[12])*algebraic[14]) rates[0] = states[1] return(rates) def computeAlgebraic(constants, states, voi): algebraic = array([[0.0] * len(voi)] * sizeAlgebraic) states = array(states) voi = array(voi) algebraic[3] = ((0.100000*(states[0]+25.0000))/(exp((states[0]+25.0000)/10.0000)-1.00000))*constants[14] algebraic[7] = 4.00000*exp(states[0]/18.0000)*constants[14] algebraic[4] = 0.0700000*exp(states[0]/20.0000)*constants[14] algebraic[8] = (1.00000/(exp((states[0]+30.0000)/10.0000)+1.00000))*constants[14] algebraic[5] = ((0.0100000*(states[0]+10.0000))/(exp((states[0]+10.0000)/10.0000)-1.00000))*constants[15] algebraic[9] = 0.125000*exp(states[0]/80.0000)*constants[15] algebraic[1] = constants[9]*(power(states[3], 3.00000))*states[2] algebraic[6] = algebraic[1]*(states[0]-constants[6]) algebraic[10] = constants[10]*(power(states[4], 4.00000)) algebraic[11] = algebraic[10]*(states[0]-constants[7]) algebraic[13] = constants[11]*(states[0]-constants[8]) algebraic[14] = algebraic[6]+algebraic[11]+algebraic[13] algebraic[0] = -constants[12]*states[1] algebraic[2] = -states[0] algebraic[12] = algebraic[10]+algebraic[1]+constants[11] algebraic[15] = (constants[12]/constants[1])*rates[1] return algebraic def solve_model(): """Solve model with ODE solver""" from scipy.integrate import ode # Initialise constants and state variables (init_states, constants) = initConsts() # Set timespan to solve over voi = linspace(0, 10, 500) # Construct ODE object to solve r = ode(computeRates) r.set_integrator('vode', method='bdf', atol=1e-06, rtol=1e-06, max_step=1) r.set_initial_value(init_states, voi[0]) r.set_f_params(constants) # Solve model states = array([[0.0] * len(voi)] * sizeStates) states[:,0] = init_states for (i,t) in enumerate(voi[1:]): if r.successful(): r.integrate(t) states[:,i+1] = r.y else: break # Compute algebraic variables algebraic = computeAlgebraic(constants, states, voi) return (voi, states, algebraic) def plot_model(voi, states, algebraic): """Plot variables against variable of integration""" import pylab (legend_states, legend_algebraic, legend_voi, legend_constants) = createLegends() pylab.figure(1) pylab.plot(voi,vstack((states,algebraic)).T) pylab.xlabel(legend_voi) pylab.legend(legend_states + legend_algebraic, loc='best') pylab.show() if __name__ == "__main__": (voi, states, algebraic) = solve_model() plot_model(voi, states, algebraic)