Generated Code
The following is python code generated by the CellML API from this CellML file. (Back to language selection)
The raw code is available.
# Size of variable arrays:
sizeAlgebraic = 88
sizeStates = 14
sizeConstants = 51
from math import *
from numpy import *
def createLegends():
legend_states = [""] * sizeStates
legend_rates = [""] * sizeStates
legend_algebraic = [""] * sizeAlgebraic
legend_voi = ""
legend_constants = [""] * sizeConstants
legend_voi = "time in component environment (millisecond)"
legend_states[0] = "Vm in component membrane (millivolt)"
legend_algebraic[76] = "I_tot in component membrane (picoampere)"
legend_constants[0] = "I_ext in component membrane (picoampere)"
legend_constants[1] = "Cm in component membrane (picofarad)"
legend_algebraic[36] = "I_Ki in component I_Ki (picoampere)"
legend_algebraic[48] = "I_CaL in component I_CaL (picoampere)"
legend_algebraic[59] = "I_VDDR in component I_VDDR (picoampere)"
legend_algebraic[64] = "I_AI in component I_AI (picoampere)"
legend_algebraic[70] = "I_NaCa in component I_NaCa (picoampere)"
legend_algebraic[72] = "I_NaK in component I_NaK (picoampere)"
legend_algebraic[73] = "I_PMCA in component I_PMCA (picoampere)"
legend_states[1] = "Na_i in component internal_ion_concentrations (millimolar)"
legend_states[2] = "K_i in component internal_ion_concentrations (millimolar)"
legend_algebraic[0] = "Ca_i in component internal_ion_concentrations (millimolar)"
legend_algebraic[2] = "Na_CF in component internal_ion_concentrations (millimolar)"
legend_algebraic[7] = "K_CF in component internal_ion_concentrations (millimolar)"
legend_algebraic[11] = "Ca_CF in component internal_ion_concentrations (millimolar)"
legend_algebraic[74] = "Inet_Na in component internal_ion_concentrations (picoampere)"
legend_algebraic[75] = "Inet_K in component internal_ion_concentrations (picoampere)"
legend_algebraic[77] = "Inet_Ca in component internal_ion_concentrations (picoampere)"
legend_states[3] = "Cai_A in component internal_ion_concentrations (millimolar)"
legend_algebraic[45] = "I_CaL_Na in component I_CaL (picoampere)"
legend_algebraic[42] = "I_CaL_K in component I_CaL (picoampere)"
legend_algebraic[39] = "I_CaL_Ca in component I_CaL (picoampere)"
legend_algebraic[63] = "I_AI_Na in component I_AI (picoampere)"
legend_algebraic[62] = "I_AI_K in component I_AI (picoampere)"
legend_algebraic[61] = "I_AI_Ca in component I_AI (picoampere)"
legend_algebraic[86] = "I_IP3R in component I_IP3R (picoampere)"
legend_algebraic[80] = "I_up in component I_up (picoampere)"
legend_algebraic[84] = "I_leak in component I_leak (picoampere)"
legend_constants[2] = "Vi in component model_parameters (micrometre3)"
legend_constants[3] = "T in component model_parameters (kelvin)"
legend_constants[4] = "R in component model_parameters (coulomb_millivolt_kelvin_millimole)"
legend_constants[5] = "F in component model_parameters (coulomb_millimole)"
legend_constants[6] = "Na_o in component model_parameters (millimolar)"
legend_constants[7] = "K_o in component model_parameters (millimolar)"
legend_constants[8] = "Ca_o in component model_parameters (millimolar)"
legend_constants[9] = "ZNa in component model_parameters (dimensionless)"
legend_constants[10] = "ZK in component model_parameters (dimensionless)"
legend_constants[11] = "ZCa in component model_parameters (dimensionless)"
legend_algebraic[20] = "fU in component I_Ki (dimensionless)"
legend_algebraic[14] = "E_K in component I_Ki (millivolt)"
legend_constants[12] = "g_Ki in component I_Ki (picoampere_per_millivolt)"
legend_algebraic[16] = "mu in component I_Ki (dimensionless)"
legend_algebraic[18] = "lambda in component I_Ki (dimensionless)"
legend_algebraic[34] = "n in component I_Ki_n_gate (dimensionless)"
legend_algebraic[22] = "n_1 in component I_Ki_n_gate (dimensionless)"
legend_algebraic[30] = "n_2 in component I_Ki_n_gate (dimensionless)"
legend_algebraic[32] = "n_3 in component I_Ki_n_gate (dimensionless)"
legend_algebraic[24] = "a in component I_Ki_n_gate (dimensionless)"
legend_algebraic[26] = "b in component I_Ki_n_gate (dimensionless)"
legend_algebraic[28] = "AB in component I_Ki_n_gate (dimensionless)"
legend_constants[13] = "PCaL in component I_CaL (picoampere_per_millimolar)"
legend_states[4] = "m in component I_CaL_m_gate (dimensionless)"
legend_states[5] = "h in component I_CaL_h_gate (dimensionless)"
legend_algebraic[1] = "alpha_m in component I_CaL_m_gate (first_order_rate_constant)"
legend_algebraic[3] = "beta_m in component I_CaL_m_gate (first_order_rate_constant)"
legend_algebraic[4] = "alpha_h in component I_CaL_h_gate (first_order_rate_constant)"
legend_algebraic[8] = "beta_h in component I_CaL_h_gate (first_order_rate_constant)"
legend_constants[14] = "PVDDR in component I_VDDR (picoampere_per_millimolar)"
legend_algebraic[53] = "m in component I_VDDR_m_gate (dimensionless)"
legend_algebraic[58] = "h in component I_VDDR_h_gate (dimensionless)"
legend_algebraic[51] = "alpha_m in component I_VDDR_m_gate (first_order_rate_constant)"
legend_algebraic[52] = "beta_m in component I_VDDR_m_gate (first_order_rate_constant)"
legend_algebraic[49] = "m2 in component I_VDDR_m_gate (first_order_rate_constant)"
legend_algebraic[50] = "m3 in component I_VDDR_m_gate (dimensionless)"
legend_algebraic[56] = "alpha_h in component I_VDDR_h_gate (first_order_rate_constant)"
legend_algebraic[57] = "beta_h in component I_VDDR_h_gate (first_order_rate_constant)"
legend_algebraic[54] = "h2 in component I_VDDR_h_gate (first_order_rate_constant)"
legend_algebraic[55] = "h3 in component I_VDDR_h_gate (dimensionless)"
legend_algebraic[60] = "po in component I_AI (dimensionless)"
legend_constants[15] = "PAI in component I_AI (picoampere_per_millimolar)"
legend_constants[16] = "Km_Ca_i in component I_AI (millimolar)"
legend_constants[17] = "PNaCa in component I_NaCa (picoampere)"
legend_constants[18] = "Km_Na_i in component I_NaCa (millimolar3)"
legend_constants[19] = "Km_Na_o in component I_NaCa (millimolar3)"
legend_constants[20] = "Km_Ca_i in component I_NaCa (millimolar)"
legend_constants[21] = "Km_Ca_o in component I_NaCa (millimolar)"
legend_constants[42] = "p_E2Na in component I_NaCa (dimensionless)"
legend_algebraic[65] = "p_E1Na in component I_NaCa (dimensionless)"
legend_algebraic[66] = "p_E1Ca in component I_NaCa (dimensionless)"
legend_constants[43] = "p_E2Ca in component I_NaCa (dimensionless)"
legend_states[6] = "y in component I_NaCa_y_gate (dimensionless)"
legend_algebraic[67] = "m in component I_NaCa_y_gate (dimensionless)"
legend_algebraic[69] = "h in component I_NaCa_y_gate (dimensionless)"
legend_algebraic[71] = "alpha_y in component I_NaCa_y_gate (first_order_rate_constant)"
legend_algebraic[68] = "beta_y in component I_NaCa_y_gate (first_order_rate_constant)"
legend_constants[22] = "pNaK in component I_NaK (picoampere)"
legend_constants[23] = "Km_Na_i in component I_NaK (millimolar)"
legend_constants[24] = "Km_K_o in component I_NaK (millimolar)"
legend_constants[25] = "pPMCA in component I_PMCA (picoampere)"
legend_constants[26] = "Km_Ca_i in component I_PMCA (millimolar)"
legend_constants[27] = "Pup in component I_up (picoampere_millisecond)"
legend_constants[28] = "Km_Ca_i in component I_up (millimolar)"
legend_constants[29] = "Km_Ca_o in component I_up (millimolar)"
legend_algebraic[83] = "p_E2 in component I_up (first_order_rate_constant)"
legend_algebraic[81] = "p_E1 in component I_up (first_order_rate_constant)"
legend_algebraic[78] = "p_E1Ca in component I_up (first_order_rate_constant)"
legend_algebraic[79] = "p_E2Ca in component I_up (first_order_rate_constant)"
legend_states[7] = "Ca_up in component calcium_concentrations_in_the_SR (millimolar)"
legend_states[8] = "y in component I_up_y_gate (dimensionless)"
legend_algebraic[85] = "alpha_y in component I_up_y_gate (first_order_rate_constant)"
legend_algebraic[87] = "beta_y in component I_up_y_gate (first_order_rate_constant)"
legend_algebraic[82] = "I_tr in component I_tr (picoampere)"
legend_constants[30] = "ptr in component I_tr (picoampere_per_millimolar)"
legend_states[9] = "Ca_rel in component calcium_concentrations_in_the_SR (millimolar)"
legend_constants[31] = "pleak in component I_leak (picoampere_per_millimolar)"
legend_constants[32] = "PIP3R in component I_IP3R (picoampere_per_millimolar)"
legend_algebraic[5] = "kco in component I_IP3R (first_order_rate_constant)"
legend_constants[33] = "koi in component I_IP3R (first_order_rate_constant)"
legend_algebraic[9] = "kic in component I_IP3R (first_order_rate_constant)"
legend_constants[34] = "kci in component I_IP3R (first_order_rate_constant)"
legend_algebraic[25] = "SC_0 in component I_IP3R (first_order_rate_constant)"
legend_algebraic[38] = "SC_3 in component I_IP3R (dimensionless)"
legend_algebraic[37] = "SC_2 in component I_IP3R (dimensionless)"
legend_algebraic[27] = "SC_1 in component I_IP3R (first_order_rate_constant)"
legend_states[10] = "m in component I_IP3R (dimensionless)"
legend_states[11] = "h in component I_IP3R (dimensionless)"
legend_algebraic[46] = "m_original in component I_IP3R (dimensionless)"
legend_algebraic[47] = "h_original in component I_IP3R (dimensionless)"
legend_algebraic[40] = "stO in component I_IP3R (dimensionless)"
legend_algebraic[41] = "stA in component I_IP3R (dimensionless)"
legend_algebraic[43] = "ndO in component I_IP3R (dimensionless)"
legend_algebraic[44] = "ndA in component I_IP3R (dimensionless)"
legend_algebraic[17] = "OrgO_N4 in component I_IP3R (first_order_rate_constant)"
legend_algebraic[19] = "OrgA_N4 in component I_IP3R (first_order_rate_constant)"
legend_algebraic[29] = "state_c1 in component I_IP3R (first_order_rate_constant)"
legend_algebraic[31] = "state_c2 in component I_IP3R (first_order_rate_constant)"
legend_algebraic[33] = "state_c3 in component I_IP3R (first_order_rate_constant)"
legend_constants[47] = "state_d1 in component I_IP3R (first_order_rate_constant)"
legend_algebraic[35] = "state_d2 in component I_IP3R (first_order_rate_constant)"
legend_constants[48] = "state_d3 in component I_IP3R (first_order_rate_constant)"
legend_algebraic[21] = "state_bf in component I_IP3R (first_order_rate_constant)"
legend_algebraic[23] = "state_sbf in component I_IP3R (per_millisecond2)"
legend_constants[41] = "kab in component I_IP3R (first_order_rate_constant)"
legend_constants[44] = "kba in component I_IP3R (first_order_rate_constant)"
legend_algebraic[12] = "kbc in component I_IP3R (first_order_rate_constant)"
legend_constants[45] = "kcb in component I_IP3R (first_order_rate_constant)"
legend_constants[46] = "kac in component I_IP3R (first_order_rate_constant)"
legend_algebraic[15] = "kca in component I_IP3R (first_order_rate_constant)"
legend_constants[35] = "dt in component I_IP3R (millisecond)"
legend_constants[36] = "dt_CellML in component I_IP3R (millisecond)"
legend_states[12] = "IP3 in component IP3_metabolism (millimolar)"
legend_constants[49] = "Vup in component model_parameters (micrometre3)"
legend_constants[50] = "Vrel in component model_parameters (micrometre3)"
legend_states[13] = "PIP2 in component IP3_metabolism (millimolar)"
legend_algebraic[6] = "PIP2_IP3 in component IP3_metabolism (first_order_rate_constant)"
legend_algebraic[10] = "IP3_PIP2 in component IP3_metabolism (first_order_rate_constant)"
legend_constants[37] = "IP3_IP4 in component IP3_metabolism (first_order_rate_constant)"
legend_constants[38] = "PIP2_IP4 in component IP3_metabolism (first_order_rate_constant)"
legend_algebraic[13] = "IP4_PIP2 in component IP3_metabolism (first_order_rate_constant)"
legend_constants[39] = "IP3total in component IP3_metabolism (millimolar)"
legend_constants[40] = "Km_Ca_i in component IP3_metabolism (millimolar)"
legend_rates[0] = "d/dt Vm in component membrane (millivolt)"
legend_rates[1] = "d/dt Na_i in component internal_ion_concentrations (millimolar)"
legend_rates[2] = "d/dt K_i in component internal_ion_concentrations (millimolar)"
legend_rates[3] = "d/dt Cai_A in component internal_ion_concentrations (millimolar)"
legend_rates[4] = "d/dt m in component I_CaL_m_gate (dimensionless)"
legend_rates[5] = "d/dt h in component I_CaL_h_gate (dimensionless)"
legend_rates[6] = "d/dt y in component I_NaCa_y_gate (dimensionless)"
legend_rates[8] = "d/dt y in component I_up_y_gate (dimensionless)"
legend_rates[10] = "d/dt m in component I_IP3R (dimensionless)"
legend_rates[11] = "d/dt h in component I_IP3R (dimensionless)"
legend_rates[7] = "d/dt Ca_up in component calcium_concentrations_in_the_SR (millimolar)"
legend_rates[9] = "d/dt Ca_rel in component calcium_concentrations_in_the_SR (millimolar)"
legend_rates[12] = "d/dt IP3 in component IP3_metabolism (millimolar)"
legend_rates[13] = "d/dt PIP2 in component IP3_metabolism (millimolar)"
return (legend_states, legend_algebraic, legend_voi, legend_constants)
def initConsts():
constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
states[0] = -68.0
constants[0] = 0.0
constants[1] = 25.0
states[1] = 5.4
states[2] = 140.0
states[3] = 0.0035
constants[2] = 715.5
constants[3] = 309.15
constants[4] = 8.314
constants[5] = 96.485
constants[6] = 140.0
constants[7] = 5.4
constants[8] = 1.8
constants[9] = 1.0
constants[10] = 1.0
constants[11] = 2.0
constants[12] = 4.32
constants[13] = 5.76
states[4] = 0.0
states[5] = 0.0
constants[14] = 7.92
constants[15] = 15.61
constants[16] = 0.003
constants[17] = 8.458
constants[18] = 670.0
constants[19] = 670000.0
constants[20] = 0.00138
constants[21] = 1.38
states[6] = 0.9459136
constants[22] = 79.0
constants[23] = 11.0
constants[24] = 0.27
constants[25] = 927.4
constants[26] = 0.008
constants[27] = 7728.0
constants[28] = 0.08
constants[29] = 0.002
states[7] = 4.6
states[8] = 0.1
constants[30] = 20.5
states[9] = 13.1
constants[31] = 0.45
constants[32] = 217.2
constants[33] = 0.02
constants[34] = 0.000849
states[10] = 0.0
states[11] = 1.0
constants[35] = 0.2
constants[36] = 0.01
states[12] = 0.0002682
states[13] = 0.002664
constants[37] = 0.004
constants[38] = 0.00
constants[39] = 0.0033
constants[40] = 0.00001
constants[41] = constants[33]
constants[42] = 1.00000/(1.00000+(1.00000+constants[8]/constants[21])/((power(constants[6], 3.00000))/constants[19]))
constants[43] = 1.00000/(1.00000+(1.00000+(power(constants[6], 3.00000))/constants[19])/(constants[8]/constants[21]))
constants[44] = 0.00000
constants[45] = constants[34]
constants[46] = 0.00000
constants[47] = -(constants[44]+constants[41]+constants[46])
constants[48] = constants[44]
constants[49] = 0.100000*constants[2]
constants[50] = 0.0100000*constants[2]
return (states, constants)
def computeRates(voi, states, constants):
rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
algebraic[1] = custom_piecewise([equal(states[0] , -30.0000), 0.00217500*((states[0]+30.0000)/(1.00000-exp((states[0]+30.0000)/-2.50000))) , True, 0.00418927])
algebraic[3] = custom_piecewise([equal(states[0] , 0.00000), 0.00157896 , True, 0.000631500*(states[0]/(exp(states[0]/2.50000)-1.00000))])
rates[4] = algebraic[1]*(1.00000-states[4])-algebraic[3]*states[4]
algebraic[4] = custom_piecewise([equal(states[0] , -34.0000), 1.00010e-05 , True, 1.77500e-06*((states[0]+34.0000)/(exp((states[0]+34.0000)/5.63300)-1.00000))])
algebraic[0] = custom_piecewise([less_equal(states[3] , 0.00000), 1.00000e-12 , True, states[3]])
algebraic[8] = 0.427000*algebraic[0]*((states[0]+64.0000)/(exp((states[0]+44.0000)/-4.16000)+1.00000))
rates[5] = algebraic[4]*(1.00000-states[5])-algebraic[8]*states[5]
algebraic[6] = 0.200000*exp((states[0]+48.5000)/18.1000)*(algebraic[0]/(algebraic[0]+constants[40]))
algebraic[10] = 0.500000*exp((states[0]+100.000)/-28.5000)
rates[12] = algebraic[6]*states[13]-(algebraic[10]+constants[37])*states[12]
algebraic[13] = 0.00350000*exp((states[0]+100.000)/-25.5000)
rates[13] = (algebraic[10]*states[12]+algebraic[13]*(constants[39]-(states[12]+states[13])))-(constants[38]+algebraic[6])*states[13]
algebraic[9] = 0.0250000*(power(states[12]/0.00100000, 3.00000))*(states[9]/1.00000)
algebraic[12] = algebraic[9]
algebraic[5] = 0.147000*(power(algebraic[0]/5.70000e-05, 3.00000))
algebraic[15] = algebraic[5]
algebraic[21] = constants[41]+constants[44]+constants[46]+algebraic[15]+algebraic[12]+constants[45]
algebraic[23] = power(algebraic[21], 2.00000)-4.00000*((algebraic[12]+constants[45]+algebraic[15])*(constants[44]+constants[41]+constants[46])-(constants[46]-algebraic[12])*(algebraic[15]+constants[44]))
algebraic[25] = (algebraic[21]+power(fabs(algebraic[23]), 0.500000))/2.00000
algebraic[29] = constants[46]-algebraic[12]
algebraic[31] = -(algebraic[12]+constants[45]+algebraic[15])
algebraic[33] = algebraic[12]
algebraic[35] = algebraic[15]-constants[44]
algebraic[37] = (algebraic[31]*constants[48]-algebraic[33]*algebraic[35])/(algebraic[29]*algebraic[35]-algebraic[31]*constants[47])
algebraic[27] = (algebraic[21]-power(fabs(algebraic[23]), 0.500000))/2.00000
algebraic[17] = algebraic[5]*states[11]-constants[33]*states[10]
algebraic[40] = ((algebraic[17]+algebraic[27]*states[10])-algebraic[37]*algebraic[27])/(algebraic[27]-algebraic[25])
algebraic[43] = ((algebraic[17]+algebraic[25]*states[10])-algebraic[37]*algebraic[25])/(algebraic[25]-algebraic[27])
algebraic[46] = algebraic[40]*exp(-algebraic[25]*constants[35])+algebraic[43]*exp(-algebraic[27]*constants[35])+algebraic[37]
rates[10] = (algebraic[46]-states[10])/constants[36]
algebraic[38] = (algebraic[29]*constants[48]-algebraic[33]*constants[47])/(algebraic[31]*constants[47]-algebraic[29]*algebraic[35])
algebraic[19] = algebraic[9]*(1.00000-(states[11]+states[10]))-(constants[34]+algebraic[5])*states[11]
algebraic[41] = ((algebraic[19]+algebraic[27]*states[11])-algebraic[38]*algebraic[27])/(algebraic[27]-algebraic[25])
algebraic[44] = ((algebraic[19]+algebraic[25]*states[11])-algebraic[38]*algebraic[25])/(algebraic[25]-algebraic[27])
algebraic[47] = algebraic[41]*exp(-algebraic[25]*constants[35])+algebraic[44]*exp(-algebraic[27]*constants[35])+algebraic[38]
rates[11] = (algebraic[47]-states[11])/constants[36]
algebraic[69] = constants[42]*exp((0.320000-1.00000)*(states[0]/((constants[4]*constants[3])/constants[5])))
algebraic[71] = 1.00000*(algebraic[69]+constants[43])
algebraic[66] = 1.00000/(1.00000+(1.00000+(power(states[1], 3.00000))/constants[18])/(algebraic[0]/constants[20]))
algebraic[65] = 1.00000/(1.00000+(1.00000+algebraic[0]/constants[20])/((power(states[1], 3.00000))/constants[18]))
algebraic[67] = algebraic[65]*exp(0.320000*(states[0]/((constants[4]*constants[3])/constants[5])))
algebraic[68] = 1.00000*(algebraic[67]+algebraic[66])
rates[6] = algebraic[71]*(1.00000-states[6])-algebraic[68]*states[6]
algebraic[70] = constants[17]*(algebraic[67]*states[6]-algebraic[69]*(1.00000-states[6]))
algebraic[72] = (constants[22]/(1.00000+power(constants[23]/states[1], 1.36000)))*((1.00000-power((states[0]+50.0000)/250.000, 2.00000))/(1.00000+constants[24]/constants[7]))
algebraic[2] = ((constants[9]*constants[5]*states[0])/(constants[4]*constants[3]))*((states[1]-constants[6]*exp((-constants[9]*constants[5]*states[0])/(constants[4]*constants[3])))/(1.00000-exp((-constants[9]*constants[5]*states[0])/(constants[4]*constants[3]))))
algebraic[45] = 5.00000e-05*constants[13]*algebraic[2]*states[4]*states[5]
algebraic[60] = 0.00170000+(1.00000-0.00170000)/(1.00000+constants[16]/algebraic[0])
algebraic[63] = 0.400000*constants[15]*algebraic[2]*algebraic[60]
algebraic[74] = algebraic[45]+algebraic[63]+3.00000*algebraic[72]+3.00000*algebraic[70]
rates[1] = -algebraic[74]/(constants[9]*constants[5]*constants[2])
algebraic[14] = ((constants[4]*constants[3])/(1.00000*constants[5]))*log(constants[7]/states[2])
algebraic[16] = 3.30000*exp((states[0]-(algebraic[14]+6.00000))/15.0000)
algebraic[18] = 26.0000*exp((states[0]-(algebraic[14]+6.00000))/52.0000)
algebraic[20] = algebraic[16]/(algebraic[16]+algebraic[18])
algebraic[22] = 1.00000/(1.00000+0.100000*exp((states[0]-algebraic[14])/15.0000)+0.0480000*exp((states[0]-algebraic[14])/7.00000))
algebraic[24] = 2.50000/(4.70000*exp((states[0]-algebraic[14])/28.7000))
algebraic[26] = 2.50000/(6.00000*exp((algebraic[14]-states[0])/25.8000))
algebraic[28] = algebraic[24]+algebraic[26]
algebraic[30] = algebraic[22]*algebraic[28]
algebraic[32] = algebraic[28]-algebraic[30]
algebraic[34] = algebraic[30]/((algebraic[30]+algebraic[32])*(power(1.00000-algebraic[20], 3.00000)))
algebraic[36] = constants[12]*(power(constants[7]/5.40000, 0.620000))*(states[0]-algebraic[14])*algebraic[34]*algebraic[20]
algebraic[7] = ((constants[10]*constants[5]*states[0])/(constants[4]*constants[3]))*((states[2]-constants[7]*exp((-constants[10]*constants[5]*states[0])/(constants[4]*constants[3])))/(1.00000-exp((-constants[10]*constants[5]*states[0])/(constants[4]*constants[3]))))
algebraic[42] = 0.00100000*constants[13]*algebraic[7]*states[4]*states[5]
algebraic[62] = 0.360000*constants[15]*algebraic[7]*algebraic[60]
algebraic[75] = (algebraic[36]+algebraic[42]+algebraic[62])-2.00000*algebraic[72]
rates[2] = -algebraic[75]/(constants[10]*constants[5]*constants[2])
algebraic[11] = ((constants[11]*constants[5]*states[0])/(constants[4]*constants[3]))*((algebraic[0]-constants[8]*exp((-constants[11]*constants[5]*states[0])/(constants[4]*constants[3])))/(1.00000-exp((-constants[11]*constants[5]*states[0])/(constants[4]*constants[3]))))
algebraic[39] = constants[13]*algebraic[11]*states[4]*states[5]
algebraic[48] = algebraic[39]+algebraic[42]+algebraic[45]
algebraic[49] = 1.00000/(1.00000+exp(-(states[0]+23.0000)/6.10000))
algebraic[50] = (4.43000-0.660000)*exp(-(power(states[0]+150.000, 2.00000))/7118.00)+0.660000
algebraic[51] = algebraic[49]/(algebraic[50]+power(10.0000, -20.0000))
algebraic[52] = (1.00000-algebraic[49])/(algebraic[50]+power(10.0000, -20.0000))
algebraic[53] = algebraic[51]/(algebraic[51]+algebraic[52])
algebraic[54] = 1.00000/(1.00000+exp(-(states[0]+75.0000)/-6.60000))
algebraic[55] = (40.8000-0.840000)*exp(-(power(states[0]+106.000, 2.00000))/2292.00)+0.840000
algebraic[56] = algebraic[54]/(algebraic[55]+power(10.0000, -20.0000))
algebraic[57] = (1.00000-algebraic[54])/(algebraic[55]+power(10.0000, -20.0000))
algebraic[58] = algebraic[56]/(algebraic[56]+algebraic[57])
algebraic[59] = constants[14]*algebraic[11]*algebraic[53]*algebraic[58]
algebraic[61] = constants[15]*algebraic[11]*algebraic[60]
algebraic[64] = algebraic[61]+algebraic[62]+algebraic[63]
algebraic[73] = constants[25]/(1.00000+power(constants[26]/algebraic[0], 2.00000))
algebraic[76] = algebraic[70]+algebraic[72]+algebraic[73]+algebraic[59]+algebraic[64]+algebraic[36]+algebraic[48]
rates[0] = -(algebraic[76]+constants[0])/constants[1]
algebraic[78] = 0.0100000/(1.00000+constants[28]/states[7])
algebraic[79] = 1.00000/(1.00000+constants[29]/algebraic[0])
algebraic[80] = 1.50000*constants[27]*(algebraic[79]*(1.00000-states[8])-algebraic[78]*states[8])
algebraic[84] = constants[31]*(states[7]-algebraic[0])
algebraic[82] = 2.80000*constants[30]*(states[7]-states[9])
rates[7] = (algebraic[80]-(algebraic[82]+algebraic[84]))/(constants[11]*constants[49]*constants[5])
algebraic[77] = (algebraic[39]+algebraic[59]+algebraic[61]+algebraic[73])-2.00000*algebraic[70]
algebraic[86] = constants[32]*(states[9]-algebraic[0])*states[10]
rates[3] = -((algebraic[77]+algebraic[80])-(algebraic[86]+algebraic[84]))/(constants[11]*constants[5]*constants[2])
algebraic[83] = 0.0100000/(1.00000+algebraic[0]/constants[29])
algebraic[85] = algebraic[79]+algebraic[83]
algebraic[81] = 1.00000/(1.00000+states[7]/constants[28])
algebraic[87] = algebraic[78]+algebraic[81]
rates[8] = algebraic[85]*(1.00000-states[8])-algebraic[87]*states[8]
rates[9] = (algebraic[82]-algebraic[86])/(constants[11]*constants[50]*constants[5])
return(rates)
def computeAlgebraic(constants, states, voi):
algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
states = array(states)
voi = array(voi)
algebraic[1] = custom_piecewise([equal(states[0] , -30.0000), 0.00217500*((states[0]+30.0000)/(1.00000-exp((states[0]+30.0000)/-2.50000))) , True, 0.00418927])
algebraic[3] = custom_piecewise([equal(states[0] , 0.00000), 0.00157896 , True, 0.000631500*(states[0]/(exp(states[0]/2.50000)-1.00000))])
algebraic[4] = custom_piecewise([equal(states[0] , -34.0000), 1.00010e-05 , True, 1.77500e-06*((states[0]+34.0000)/(exp((states[0]+34.0000)/5.63300)-1.00000))])
algebraic[0] = custom_piecewise([less_equal(states[3] , 0.00000), 1.00000e-12 , True, states[3]])
algebraic[8] = 0.427000*algebraic[0]*((states[0]+64.0000)/(exp((states[0]+44.0000)/-4.16000)+1.00000))
algebraic[6] = 0.200000*exp((states[0]+48.5000)/18.1000)*(algebraic[0]/(algebraic[0]+constants[40]))
algebraic[10] = 0.500000*exp((states[0]+100.000)/-28.5000)
algebraic[13] = 0.00350000*exp((states[0]+100.000)/-25.5000)
algebraic[9] = 0.0250000*(power(states[12]/0.00100000, 3.00000))*(states[9]/1.00000)
algebraic[12] = algebraic[9]
algebraic[5] = 0.147000*(power(algebraic[0]/5.70000e-05, 3.00000))
algebraic[15] = algebraic[5]
algebraic[21] = constants[41]+constants[44]+constants[46]+algebraic[15]+algebraic[12]+constants[45]
algebraic[23] = power(algebraic[21], 2.00000)-4.00000*((algebraic[12]+constants[45]+algebraic[15])*(constants[44]+constants[41]+constants[46])-(constants[46]-algebraic[12])*(algebraic[15]+constants[44]))
algebraic[25] = (algebraic[21]+power(fabs(algebraic[23]), 0.500000))/2.00000
algebraic[29] = constants[46]-algebraic[12]
algebraic[31] = -(algebraic[12]+constants[45]+algebraic[15])
algebraic[33] = algebraic[12]
algebraic[35] = algebraic[15]-constants[44]
algebraic[37] = (algebraic[31]*constants[48]-algebraic[33]*algebraic[35])/(algebraic[29]*algebraic[35]-algebraic[31]*constants[47])
algebraic[27] = (algebraic[21]-power(fabs(algebraic[23]), 0.500000))/2.00000
algebraic[17] = algebraic[5]*states[11]-constants[33]*states[10]
algebraic[40] = ((algebraic[17]+algebraic[27]*states[10])-algebraic[37]*algebraic[27])/(algebraic[27]-algebraic[25])
algebraic[43] = ((algebraic[17]+algebraic[25]*states[10])-algebraic[37]*algebraic[25])/(algebraic[25]-algebraic[27])
algebraic[46] = algebraic[40]*exp(-algebraic[25]*constants[35])+algebraic[43]*exp(-algebraic[27]*constants[35])+algebraic[37]
algebraic[38] = (algebraic[29]*constants[48]-algebraic[33]*constants[47])/(algebraic[31]*constants[47]-algebraic[29]*algebraic[35])
algebraic[19] = algebraic[9]*(1.00000-(states[11]+states[10]))-(constants[34]+algebraic[5])*states[11]
algebraic[41] = ((algebraic[19]+algebraic[27]*states[11])-algebraic[38]*algebraic[27])/(algebraic[27]-algebraic[25])
algebraic[44] = ((algebraic[19]+algebraic[25]*states[11])-algebraic[38]*algebraic[25])/(algebraic[25]-algebraic[27])
algebraic[47] = algebraic[41]*exp(-algebraic[25]*constants[35])+algebraic[44]*exp(-algebraic[27]*constants[35])+algebraic[38]
algebraic[69] = constants[42]*exp((0.320000-1.00000)*(states[0]/((constants[4]*constants[3])/constants[5])))
algebraic[71] = 1.00000*(algebraic[69]+constants[43])
algebraic[66] = 1.00000/(1.00000+(1.00000+(power(states[1], 3.00000))/constants[18])/(algebraic[0]/constants[20]))
algebraic[65] = 1.00000/(1.00000+(1.00000+algebraic[0]/constants[20])/((power(states[1], 3.00000))/constants[18]))
algebraic[67] = algebraic[65]*exp(0.320000*(states[0]/((constants[4]*constants[3])/constants[5])))
algebraic[68] = 1.00000*(algebraic[67]+algebraic[66])
algebraic[70] = constants[17]*(algebraic[67]*states[6]-algebraic[69]*(1.00000-states[6]))
algebraic[72] = (constants[22]/(1.00000+power(constants[23]/states[1], 1.36000)))*((1.00000-power((states[0]+50.0000)/250.000, 2.00000))/(1.00000+constants[24]/constants[7]))
algebraic[2] = ((constants[9]*constants[5]*states[0])/(constants[4]*constants[3]))*((states[1]-constants[6]*exp((-constants[9]*constants[5]*states[0])/(constants[4]*constants[3])))/(1.00000-exp((-constants[9]*constants[5]*states[0])/(constants[4]*constants[3]))))
algebraic[45] = 5.00000e-05*constants[13]*algebraic[2]*states[4]*states[5]
algebraic[60] = 0.00170000+(1.00000-0.00170000)/(1.00000+constants[16]/algebraic[0])
algebraic[63] = 0.400000*constants[15]*algebraic[2]*algebraic[60]
algebraic[74] = algebraic[45]+algebraic[63]+3.00000*algebraic[72]+3.00000*algebraic[70]
algebraic[14] = ((constants[4]*constants[3])/(1.00000*constants[5]))*log(constants[7]/states[2])
algebraic[16] = 3.30000*exp((states[0]-(algebraic[14]+6.00000))/15.0000)
algebraic[18] = 26.0000*exp((states[0]-(algebraic[14]+6.00000))/52.0000)
algebraic[20] = algebraic[16]/(algebraic[16]+algebraic[18])
algebraic[22] = 1.00000/(1.00000+0.100000*exp((states[0]-algebraic[14])/15.0000)+0.0480000*exp((states[0]-algebraic[14])/7.00000))
algebraic[24] = 2.50000/(4.70000*exp((states[0]-algebraic[14])/28.7000))
algebraic[26] = 2.50000/(6.00000*exp((algebraic[14]-states[0])/25.8000))
algebraic[28] = algebraic[24]+algebraic[26]
algebraic[30] = algebraic[22]*algebraic[28]
algebraic[32] = algebraic[28]-algebraic[30]
algebraic[34] = algebraic[30]/((algebraic[30]+algebraic[32])*(power(1.00000-algebraic[20], 3.00000)))
algebraic[36] = constants[12]*(power(constants[7]/5.40000, 0.620000))*(states[0]-algebraic[14])*algebraic[34]*algebraic[20]
algebraic[7] = ((constants[10]*constants[5]*states[0])/(constants[4]*constants[3]))*((states[2]-constants[7]*exp((-constants[10]*constants[5]*states[0])/(constants[4]*constants[3])))/(1.00000-exp((-constants[10]*constants[5]*states[0])/(constants[4]*constants[3]))))
algebraic[42] = 0.00100000*constants[13]*algebraic[7]*states[4]*states[5]
algebraic[62] = 0.360000*constants[15]*algebraic[7]*algebraic[60]
algebraic[75] = (algebraic[36]+algebraic[42]+algebraic[62])-2.00000*algebraic[72]
algebraic[11] = ((constants[11]*constants[5]*states[0])/(constants[4]*constants[3]))*((algebraic[0]-constants[8]*exp((-constants[11]*constants[5]*states[0])/(constants[4]*constants[3])))/(1.00000-exp((-constants[11]*constants[5]*states[0])/(constants[4]*constants[3]))))
algebraic[39] = constants[13]*algebraic[11]*states[4]*states[5]
algebraic[48] = algebraic[39]+algebraic[42]+algebraic[45]
algebraic[49] = 1.00000/(1.00000+exp(-(states[0]+23.0000)/6.10000))
algebraic[50] = (4.43000-0.660000)*exp(-(power(states[0]+150.000, 2.00000))/7118.00)+0.660000
algebraic[51] = algebraic[49]/(algebraic[50]+power(10.0000, -20.0000))
algebraic[52] = (1.00000-algebraic[49])/(algebraic[50]+power(10.0000, -20.0000))
algebraic[53] = algebraic[51]/(algebraic[51]+algebraic[52])
algebraic[54] = 1.00000/(1.00000+exp(-(states[0]+75.0000)/-6.60000))
algebraic[55] = (40.8000-0.840000)*exp(-(power(states[0]+106.000, 2.00000))/2292.00)+0.840000
algebraic[56] = algebraic[54]/(algebraic[55]+power(10.0000, -20.0000))
algebraic[57] = (1.00000-algebraic[54])/(algebraic[55]+power(10.0000, -20.0000))
algebraic[58] = algebraic[56]/(algebraic[56]+algebraic[57])
algebraic[59] = constants[14]*algebraic[11]*algebraic[53]*algebraic[58]
algebraic[61] = constants[15]*algebraic[11]*algebraic[60]
algebraic[64] = algebraic[61]+algebraic[62]+algebraic[63]
algebraic[73] = constants[25]/(1.00000+power(constants[26]/algebraic[0], 2.00000))
algebraic[76] = algebraic[70]+algebraic[72]+algebraic[73]+algebraic[59]+algebraic[64]+algebraic[36]+algebraic[48]
algebraic[78] = 0.0100000/(1.00000+constants[28]/states[7])
algebraic[79] = 1.00000/(1.00000+constants[29]/algebraic[0])
algebraic[80] = 1.50000*constants[27]*(algebraic[79]*(1.00000-states[8])-algebraic[78]*states[8])
algebraic[84] = constants[31]*(states[7]-algebraic[0])
algebraic[82] = 2.80000*constants[30]*(states[7]-states[9])
algebraic[77] = (algebraic[39]+algebraic[59]+algebraic[61]+algebraic[73])-2.00000*algebraic[70]
algebraic[86] = constants[32]*(states[9]-algebraic[0])*states[10]
algebraic[83] = 0.0100000/(1.00000+algebraic[0]/constants[29])
algebraic[85] = algebraic[79]+algebraic[83]
algebraic[81] = 1.00000/(1.00000+states[7]/constants[28])
algebraic[87] = algebraic[78]+algebraic[81]
return algebraic
def custom_piecewise(cases):
"""Compute result of a piecewise function"""
return select(cases[0::2],cases[1::2])
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)