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 = 66
sizeStates = 35
sizeConstants = 76
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] = "V in component membrane (millivolt)"
legend_constants[0] = "Cm in component membrane (microF_per_cm2)"
legend_constants[1] = "Vmyo in component membrane (microlitre)"
legend_constants[2] = "VJSR in component membrane (microlitre)"
legend_constants[3] = "VNSR in component membrane (microlitre)"
legend_constants[4] = "Vss in component membrane (microlitre)"
legend_constants[5] = "Acap in component membrane (cm2)"
legend_constants[6] = "Ko in component membrane (micromolar)"
legend_constants[7] = "Nao in component membrane (micromolar)"
legend_constants[8] = "Cao in component membrane (micromolar)"
legend_constants[9] = "Clo in component membrane (micromolar)"
legend_constants[10] = "Cli in component membrane (micromolar)"
legend_constants[11] = "R in component membrane (joule_per_mole_kelvin)"
legend_constants[12] = "T in component membrane (kelvin)"
legend_constants[13] = "F in component membrane (coulomb_per_millimole)"
legend_algebraic[0] = "i_stim in component membrane (picoA_per_picoF)"
legend_constants[14] = "stim_start in component membrane (millisecond)"
legend_constants[15] = "stim_end in component membrane (millisecond)"
legend_constants[16] = "stim_period in component membrane (millisecond)"
legend_constants[17] = "stim_duration in component membrane (millisecond)"
legend_constants[18] = "stim_amplitude in component membrane (picoA_per_picoF)"
legend_algebraic[43] = "i_CaL in component L_type_calcium_current (picoA_per_picoF)"
legend_algebraic[45] = "i_CaT in component T_type_calcium_current (picoA_per_picoF)"
legend_algebraic[47] = "i_pCa in component calcium_pump_current (picoA_per_picoF)"
legend_algebraic[49] = "i_NaCa in component sodium_calcium_exchange_current (picoA_per_picoF)"
legend_algebraic[51] = "i_Cab in component calcium_background_current (picoA_per_picoF)"
legend_algebraic[53] = "i_Na in component fast_sodium_current (picoA_per_picoF)"
legend_algebraic[54] = "i_Nab in component sodium_background_current (picoA_per_picoF)"
legend_algebraic[63] = "i_NaK in component sodium_potassium_pump_current (picoA_per_picoF)"
legend_algebraic[56] = "i_Kto_f in component fast_transient_outward_potassium_current (picoA_per_picoF)"
legend_algebraic[59] = "i_K1 in component time_independent_potassium_current (picoA_per_picoF)"
legend_algebraic[60] = "i_Ks in component slow_delayed_rectifier_potassium_current (picoA_per_picoF)"
legend_algebraic[57] = "i_Kur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (picoA_per_picoF)"
legend_algebraic[58] = "i_Kss in component non_inactivating_steady_state_potassium_current (picoA_per_picoF)"
legend_algebraic[65] = "i_ClCa in component calcium_activated_chloride_current (picoA_per_picoF)"
legend_algebraic[61] = "i_Kr in component rapid_delayed_rectifier_potassium_current (picoA_per_picoF)"
legend_states[1] = "Cai in component calcium_concentration (micromolar)"
legend_states[2] = "Cass in component calcium_concentration (micromolar)"
legend_states[3] = "CaJSR in component calcium_concentration (micromolar)"
legend_states[4] = "CaNSR in component calcium_concentration (micromolar)"
legend_algebraic[15] = "Bi in component calcium_concentration (dimensionless)"
legend_algebraic[30] = "Bss in component calcium_concentration (dimensionless)"
legend_algebraic[35] = "BJSR in component calcium_concentration (dimensionless)"
legend_constants[19] = "CMDN_tot in component calcium_concentration (micromolar)"
legend_constants[20] = "CSQN_tot in component calcium_concentration (micromolar)"
legend_constants[21] = "Km_CMDN in component calcium_concentration (micromolar)"
legend_constants[22] = "Km_CSQN in component calcium_concentration (micromolar)"
legend_algebraic[44] = "J_leak in component calcium_fluxes (micromolar_per_millisecond)"
legend_algebraic[39] = "J_rel in component calcium_fluxes (micromolar_per_millisecond)"
legend_algebraic[46] = "J_up in component calcium_fluxes (micromolar_per_millisecond)"
legend_algebraic[41] = "J_tr in component calcium_fluxes (micromolar_per_millisecond)"
legend_algebraic[48] = "J_trpn in component calcium_fluxes (micromolar_per_millisecond)"
legend_algebraic[42] = "J_xfer in component calcium_fluxes (micromolar_per_millisecond)"
legend_constants[23] = "k_plus_htrpn in component calcium_fluxes (per_micromolar_millisecond)"
legend_constants[24] = "k_minus_htrpn in component calcium_fluxes (per_millisecond)"
legend_constants[25] = "k_plus_ltrpn in component calcium_fluxes (per_micromolar_millisecond)"
legend_constants[26] = "k_minus_ltrpn in component calcium_fluxes (per_millisecond)"
legend_states[5] = "P_RyR in component calcium_fluxes (dimensionless)"
legend_constants[27] = "v1 in component calcium_fluxes (per_millisecond)"
legend_constants[28] = "tau_tr in component calcium_fluxes (millisecond)"
legend_constants[29] = "v2 in component calcium_fluxes (per_millisecond)"
legend_constants[30] = "tau_xfer in component calcium_fluxes (millisecond)"
legend_constants[31] = "v3 in component calcium_fluxes (micromolar_per_millisecond)"
legend_constants[32] = "Km_up in component calcium_fluxes (micromolar)"
legend_constants[33] = "LTRPN_tot in component calcium_buffering (micromolar)"
legend_constants[34] = "HTRPN_tot in component calcium_buffering (micromolar)"
legend_states[6] = "LTRPN_Ca in component calcium_buffering (micromolar)"
legend_states[7] = "HTRPN_Ca in component calcium_buffering (micromolar)"
legend_constants[35] = "i_CaL_max in component L_type_calcium_current (picoA_per_picoF)"
legend_states[8] = "P_O1 in component ryanodine_receptors (dimensionless)"
legend_states[9] = "P_O2 in component ryanodine_receptors (dimensionless)"
legend_algebraic[1] = "P_C1 in component ryanodine_receptors (dimensionless)"
legend_states[10] = "P_C2 in component ryanodine_receptors (dimensionless)"
legend_constants[36] = "k_plus_a in component ryanodine_receptors (micromolar4_per_millisecond)"
legend_constants[37] = "k_minus_a in component ryanodine_receptors (per_millisecond)"
legend_constants[38] = "k_plus_b in component ryanodine_receptors (micromolar3_per_millisecond)"
legend_constants[39] = "k_minus_b in component ryanodine_receptors (per_millisecond)"
legend_constants[40] = "k_plus_c in component ryanodine_receptors (per_millisecond)"
legend_constants[41] = "k_minus_c in component ryanodine_receptors (per_millisecond)"
legend_constants[42] = "m in component ryanodine_receptors (dimensionless)"
legend_constants[43] = "n in component ryanodine_receptors (dimensionless)"
legend_constants[44] = "E_CaL in component L_type_calcium_current (millivolt)"
legend_constants[45] = "g_CaL in component L_type_calcium_current (milliS_per_microF)"
legend_states[11] = "O in component L_type_calcium_current (dimensionless)"
legend_algebraic[2] = "C1 in component L_type_calcium_current (dimensionless)"
legend_states[12] = "C2 in component L_type_calcium_current (dimensionless)"
legend_states[13] = "C3 in component L_type_calcium_current (dimensionless)"
legend_states[14] = "C4 in component L_type_calcium_current (dimensionless)"
legend_states[15] = "I1 in component L_type_calcium_current (dimensionless)"
legend_states[16] = "I2 in component L_type_calcium_current (dimensionless)"
legend_states[17] = "I3 in component L_type_calcium_current (dimensionless)"
legend_algebraic[16] = "alpha in component L_type_calcium_current (per_millisecond)"
legend_algebraic[31] = "beta in component L_type_calcium_current (per_millisecond)"
legend_algebraic[36] = "gamma in component L_type_calcium_current (per_millisecond)"
legend_algebraic[40] = "Kpcf in component L_type_calcium_current (per_millisecond)"
legend_constants[46] = "Kpcb in component L_type_calcium_current (per_millisecond)"
legend_constants[47] = "Kpc_max in component L_type_calcium_current (per_millisecond)"
legend_constants[48] = "Kpc_half in component L_type_calcium_current (micromolar)"
legend_constants[49] = "E_CaT in component T_type_calcium_current (millivolt)"
legend_constants[50] = "g_CaT in component T_type_calcium_current (milliS_per_microF)"
legend_states[18] = "b in component T_type_calcium_current_b_gate (dimensionless)"
legend_states[19] = "g in component T_type_calcium_current_g_gate (dimensionless)"
legend_algebraic[3] = "b_infinity in component T_type_calcium_current_b_gate (dimensionless)"
legend_algebraic[17] = "tau_b in component T_type_calcium_current_b_gate (millisecond)"
legend_algebraic[4] = "g_infinity in component T_type_calcium_current_g_gate (dimensionless)"
legend_algebraic[18] = "tau_g in component T_type_calcium_current_g_gate (millisecond)"
legend_constants[51] = "i_pCa_max in component calcium_pump_current (picoA_per_picoF)"
legend_constants[52] = "Km_pCa in component calcium_pump_current (micromolar)"
legend_constants[53] = "k_NaCa in component sodium_calcium_exchange_current (picoA_per_picoF)"
legend_constants[54] = "K_mNa in component sodium_calcium_exchange_current (micromolar)"
legend_constants[55] = "K_mCa in component sodium_calcium_exchange_current (micromolar)"
legend_constants[56] = "k_sat in component sodium_calcium_exchange_current (dimensionless)"
legend_constants[57] = "eta in component sodium_calcium_exchange_current (dimensionless)"
legend_states[20] = "Nai in component sodium_concentration (micromolar)"
legend_constants[58] = "g_Cab in component calcium_background_current (milliS_per_microF)"
legend_algebraic[50] = "E_CaN in component calcium_background_current (millivolt)"
legend_algebraic[52] = "E_Na in component fast_sodium_current (millivolt)"
legend_constants[59] = "g_Na in component fast_sodium_current (milliS_per_microF)"
legend_states[21] = "m in component fast_sodium_current_m_gate (dimensionless)"
legend_states[22] = "h in component fast_sodium_current_h_gate (dimensionless)"
legend_states[23] = "j in component fast_sodium_current_j_gate (dimensionless)"
legend_algebraic[5] = "alpha_m in component fast_sodium_current_m_gate (per_millisecond)"
legend_algebraic[19] = "beta_m in component fast_sodium_current_m_gate (per_millisecond)"
legend_algebraic[6] = "alpha_h in component fast_sodium_current_h_gate (per_millisecond)"
legend_algebraic[20] = "beta_h in component fast_sodium_current_h_gate (per_millisecond)"
legend_algebraic[7] = "alpha_j in component fast_sodium_current_j_gate (per_millisecond)"
legend_algebraic[21] = "beta_j in component fast_sodium_current_j_gate (per_millisecond)"
legend_constants[60] = "g_Nab in component sodium_background_current (milliS_per_microF)"
legend_states[24] = "Ki in component potassium_concentration (micromolar)"
legend_algebraic[55] = "E_K in component fast_transient_outward_potassium_current (millivolt)"
legend_constants[61] = "g_Kto_f in component fast_transient_outward_potassium_current (milliS_per_microF)"
legend_states[25] = "ato_f in component fast_transient_outward_potassium_current (dimensionless)"
legend_states[26] = "ito_f in component fast_transient_outward_potassium_current (dimensionless)"
legend_algebraic[37] = "ito_f_infinity in component fast_transient_outward_potassium_current (dimensionless)"
legend_algebraic[8] = "tau_ito_f in component fast_transient_outward_potassium_current (millisecond)"
legend_algebraic[9] = "alpha_a in component fast_transient_outward_potassium_current (per_millisecond)"
legend_algebraic[23] = "beta_a in component fast_transient_outward_potassium_current (per_millisecond)"
legend_algebraic[22] = "alpha_i in component fast_transient_outward_potassium_current (per_millisecond)"
legend_algebraic[32] = "beta_i in component fast_transient_outward_potassium_current (per_millisecond)"
legend_algebraic[10] = "ass in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)"
legend_constants[62] = "g_Kur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (milliS_per_microF)"
legend_states[27] = "aur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)"
legend_states[28] = "iur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)"
legend_algebraic[24] = "tau_aur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (millisecond)"
legend_algebraic[25] = "tau_iur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (millisecond)"
legend_algebraic[11] = "iss in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)"
legend_constants[63] = "g_Kss in component non_inactivating_steady_state_potassium_current (milliS_per_microF)"
legend_states[29] = "aKss in component non_inactivating_steady_state_potassium_current (dimensionless)"
legend_algebraic[26] = "tau_Kss in component non_inactivating_steady_state_potassium_current (millisecond)"
legend_constants[64] = "g_K1 in component time_independent_potassium_current (milliS_per_microF)"
legend_constants[65] = "g_Ks in component slow_delayed_rectifier_potassium_current (milliS_per_microF)"
legend_states[30] = "nKs in component slow_delayed_rectifier_potassium_current (dimensionless)"
legend_algebraic[12] = "alpha_n in component slow_delayed_rectifier_potassium_current (per_millisecond)"
legend_algebraic[27] = "beta_n in component slow_delayed_rectifier_potassium_current (per_millisecond)"
legend_constants[66] = "g_Kr in component rapid_delayed_rectifier_potassium_current (milliS_per_microF)"
legend_states[31] = "O_K in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_states[32] = "C_K1 in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_states[33] = "C_K2 in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_algebraic[13] = "C_K0 in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_states[34] = "I_K in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_algebraic[28] = "alpha_a0 in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_algebraic[33] = "beta_a0 in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_constants[67] = "kb in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_constants[68] = "kf in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_algebraic[14] = "alpha_a1 in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_algebraic[29] = "beta_a1 in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_algebraic[34] = "alpha_i in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_algebraic[38] = "beta_i in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_constants[69] = "i_NaK_max in component sodium_potassium_pump_current (picoA_per_picoF)"
legend_constants[70] = "Km_Nai in component sodium_potassium_pump_current (micromolar)"
legend_constants[71] = "Km_Ko in component sodium_potassium_pump_current (micromolar)"
legend_algebraic[62] = "f_NaK in component sodium_potassium_pump_current (dimensionless)"
legend_constants[75] = "sigma in component sodium_potassium_pump_current (dimensionless)"
legend_constants[72] = "P_ClCa in component calcium_activated_chloride_current (cm_per_second)"
legend_algebraic[64] = "f_ClCa in component calcium_activated_chloride_current (dimensionless)"
legend_constants[73] = "Km_Cl in component calcium_activated_chloride_current (micromolar)"
legend_constants[74] = "n_ClCa in component calcium_activated_chloride_current (dimensionless)"
legend_rates[0] = "d/dt V in component membrane (millivolt)"
legend_rates[1] = "d/dt Cai in component calcium_concentration (micromolar)"
legend_rates[2] = "d/dt Cass in component calcium_concentration (micromolar)"
legend_rates[3] = "d/dt CaJSR in component calcium_concentration (micromolar)"
legend_rates[4] = "d/dt CaNSR in component calcium_concentration (micromolar)"
legend_rates[5] = "d/dt P_RyR in component calcium_fluxes (dimensionless)"
legend_rates[6] = "d/dt LTRPN_Ca in component calcium_buffering (micromolar)"
legend_rates[7] = "d/dt HTRPN_Ca in component calcium_buffering (micromolar)"
legend_rates[8] = "d/dt P_O1 in component ryanodine_receptors (dimensionless)"
legend_rates[9] = "d/dt P_O2 in component ryanodine_receptors (dimensionless)"
legend_rates[10] = "d/dt P_C2 in component ryanodine_receptors (dimensionless)"
legend_rates[11] = "d/dt O in component L_type_calcium_current (dimensionless)"
legend_rates[12] = "d/dt C2 in component L_type_calcium_current (dimensionless)"
legend_rates[13] = "d/dt C3 in component L_type_calcium_current (dimensionless)"
legend_rates[14] = "d/dt C4 in component L_type_calcium_current (dimensionless)"
legend_rates[15] = "d/dt I1 in component L_type_calcium_current (dimensionless)"
legend_rates[16] = "d/dt I2 in component L_type_calcium_current (dimensionless)"
legend_rates[17] = "d/dt I3 in component L_type_calcium_current (dimensionless)"
legend_rates[18] = "d/dt b in component T_type_calcium_current_b_gate (dimensionless)"
legend_rates[19] = "d/dt g in component T_type_calcium_current_g_gate (dimensionless)"
legend_rates[20] = "d/dt Nai in component sodium_concentration (micromolar)"
legend_rates[21] = "d/dt m in component fast_sodium_current_m_gate (dimensionless)"
legend_rates[22] = "d/dt h in component fast_sodium_current_h_gate (dimensionless)"
legend_rates[23] = "d/dt j in component fast_sodium_current_j_gate (dimensionless)"
legend_rates[24] = "d/dt Ki in component potassium_concentration (micromolar)"
legend_rates[25] = "d/dt ato_f in component fast_transient_outward_potassium_current (dimensionless)"
legend_rates[26] = "d/dt ito_f in component fast_transient_outward_potassium_current (dimensionless)"
legend_rates[27] = "d/dt aur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)"
legend_rates[28] = "d/dt iur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)"
legend_rates[29] = "d/dt aKss in component non_inactivating_steady_state_potassium_current (dimensionless)"
legend_rates[30] = "d/dt nKs in component slow_delayed_rectifier_potassium_current (dimensionless)"
legend_rates[33] = "d/dt C_K2 in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_rates[32] = "d/dt C_K1 in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_rates[31] = "d/dt O_K in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_rates[34] = "d/dt I_K in component rapid_delayed_rectifier_potassium_current (dimensionless)"
return (legend_states, legend_algebraic, legend_voi, legend_constants)
def initConsts():
constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
states[0] = -80.6475
constants[0] = 1
constants[1] = 2.2826e-6
constants[2] = 0.12e-8
constants[3] = 2.098e-7
constants[4] = 3.0734e-8
constants[5] = 1.5410e-5
constants[6] = 5400
constants[7] = 140000
constants[8] = 1000
constants[9] = 132000
constants[10] = 30000
constants[11] = 8.314
constants[12] = 298
constants[13] = 96.5
constants[14] = 20
constants[15] = 500000
constants[16] = 2000.0
constants[17] = 0.5
constants[18] = -80
states[1] = 0.2049
states[2] = 0.1867
states[3] = 705.5417
states[4] = 703.3063
constants[19] = 25
constants[20] = 15000
constants[21] = 0.238
constants[22] = 800
constants[23] = 0.00237
constants[24] = 3.2e-5
constants[25] = 0.0327
constants[26] = 0.0196
states[5] = 2.8603e-14
constants[27] = 0.45
constants[28] = 20
constants[29] = 2.088e-5
constants[30] = 8
constants[31] = 0.09
constants[32] = 0.5
constants[33] = 35
constants[34] = 70
states[6] = 8.9220
states[7] = 66.0407
constants[35] = 7
states[8] = 1.1e-3
states[9] = 2.9799e-8
states[10] = 0.1003
constants[36] = 0.00608
constants[37] = 0.07125
constants[38] = 0.00405
constants[39] = 0.965
constants[40] = 0.009
constants[41] = 0.0008
constants[42] = 3
constants[43] = 4
constants[44] = 63
constants[45] = 0.19019
states[11] = 4.4776e-12
states[12] = 1.6996e-4
states[13] = 1.1754e-8
states[14] = 6.1239e-10
states[15] = 1.5450e-9
states[16] = 6.4226e-8
states[17] = 6.5201e-7
constants[46] = 0.0005
constants[47] = 0.23324
constants[48] = 20
constants[49] = 50
constants[50] = 0.055
states[18] = 0.0047
states[19] = 0.9020
constants[51] = 0.2
constants[52] = 0.5
constants[53] = 907.68
constants[54] = 87500
constants[55] = 1380
constants[56] = 0.1
constants[57] = 0.35
states[20] = 2.1747e4
constants[58] = 0.00025
constants[59] = 10
states[21] = 0.0032
states[22] = 0.9606
states[23] = 0.9744
constants[60] = 0.0026
states[24] = 1.3645e5
constants[61] = 0.1017
states[25] = 0.0032
states[26] = 0.7530
constants[62] = 0.0048
states[27] = 5.2471e-4
states[28] = 0.9866
constants[63] = 0.015
states[29] = 0.8500
constants[64] = 0.235
constants[65] = 0.046
states[30] = 0.0026
constants[66] = 01.17
states[31] = 2.1e-3
states[32] = 1.1e-3
states[33] = 8.3641e-4
states[34] = 4.3522e-4
constants[67] = 0.036778
constants[68] = 0.023761
constants[69] = 0.88
constants[70] = 21000
constants[71] = 1500
constants[72] = 2.74e-7
constants[73] = 4
constants[74] = 3
constants[75] = (1.00000/7.00000)*(exp(constants[7]/67300.0)-1.00000)
return (states, constants)
def computeRates(voi, states, constants):
rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
rates[6] = constants[25]*states[1]*(constants[33]-states[6])-constants[26]*states[6]
rates[7] = constants[23]*states[1]*(constants[34]-states[7])-constants[24]*states[7]
rates[9] = constants[38]*(power(states[2], constants[42]))*states[8]-constants[39]*states[9]
rates[10] = constants[40]*states[8]-constants[41]*states[10]
algebraic[1] = 1.00000-(states[10]+states[8]+states[9])
rates[8] = (constants[36]*(power(states[2], constants[43]))*algebraic[1]+constants[39]*states[9]+constants[41]*states[10])-(constants[37]*states[8]+constants[38]*(power(states[2], constants[42]))*states[8]+constants[40]*states[8])
algebraic[3] = 1.00000/(1.00000+exp(-(states[0]+48.0000)/6.10000))
algebraic[17] = 0.100000+5.40000/(1.00000+exp((states[0]+100.000)/6.60000))
rates[18] = (algebraic[3]-states[18])/algebraic[17]
algebraic[4] = 1.00000/(1.00000+exp((states[0]+66.0000)/6.60000))
algebraic[18] = 8.00000+32.0000/(1.00000+exp((states[0]+65.0000)/5.00000))
rates[19] = (algebraic[4]-states[19])/algebraic[18]
algebraic[5] = (0.320000*(states[0]+47.1300))/(1.00000-exp(-0.100000*(states[0]+47.1300)))
algebraic[19] = 0.0800000*exp(-states[0]/11.0000)
rates[21] = algebraic[5]-(algebraic[5]+algebraic[19])*states[21]
algebraic[6] = custom_piecewise([less(states[0] , -40.0000), 0.135000*exp((80.0000+states[0])/-6.80000) , True, 0.00000])
algebraic[20] = custom_piecewise([less(states[0] , -40.0000), 3.56000*exp(0.0790000*states[0])+310000.*exp(0.350000*states[0]) , True, 1.00000/(0.130000*(1.00000+exp((states[0]+10.6600)/-11.1000)))])
rates[22] = algebraic[6]-(algebraic[6]+algebraic[20])*states[22]
algebraic[7] = custom_piecewise([less(states[0] , -40.0000), (-(127140.*exp(0.244400*states[0])+3.47400e-05*exp(-0.0439100*states[0]))*(states[0]+37.7800))/(1.00000+exp(0.311000*(states[0]+79.2300))) , True, 0.00000])
algebraic[21] = custom_piecewise([less(states[0] , -40.0000), (0.121200*exp(-0.0105200*states[0]))/(1.00000+exp(-0.137800*(states[0]+40.1400))) , True, (0.300000*exp(-2.53500e-07*states[0]))/(1.00000+exp(-0.100000*(states[0]+32.0000)))])
rates[23] = algebraic[7]-(algebraic[7]+algebraic[21])*states[23]
algebraic[9] = 0.180640*exp(0.0357700*(states[0]+30.0000))
algebraic[23] = 0.395600*exp(-0.0623700*(states[0]+30.0000))
rates[25] = algebraic[9]*(1.00000-states[25])-algebraic[23]*states[25]
algebraic[10] = 1.00000/(1.00000+exp(-(states[0]+22.5000)/7.70000))
algebraic[24] = 0.493000*exp(-0.0629000*states[0])+2.05800
rates[27] = (algebraic[10]-states[27])/algebraic[24]
algebraic[25] = 1200.00-170.000/(1.00000+exp((states[0]+45.2000)/5.70000))
algebraic[11] = 1.00000/(1.00000+exp((states[0]+45.2000)/5.70000))
rates[28] = (algebraic[11]-states[28])/algebraic[25]
algebraic[26] = 39.3000*exp(-0.0862000*states[0])+13.1700
rates[29] = (algebraic[10]-states[29])/algebraic[26]
algebraic[12] = (4.81333e-06*(states[0]+26.5000))/(1.00000-exp(-0.128000*(states[0]+26.5000)))
algebraic[27] = 9.53333e-05*exp(-0.0380000*(states[0]+26.5000))
rates[30] = algebraic[12]*(1.00000-states[30])-algebraic[27]*states[30]
algebraic[14] = 0.0137330*exp(0.0381980*states[0])
algebraic[29] = 6.89000e-05*exp(-0.0417800*states[0])
rates[33] = (constants[68]*states[32]+algebraic[29]*states[31])-(constants[67]*states[33]+algebraic[14]*states[33])
algebraic[2] = 1.00000-(states[11]+states[12]+states[12]+states[13]+states[14]+states[15]+states[16]+states[17])
algebraic[16] = (0.400000*exp((states[0]+12.0000)/10.0000)*((1.00000+0.700000*exp(-(power(states[0]+40.0000, 2.00000))/10.0000))-0.750000*exp(-(power(states[0]+20.0000, 2.00000))/400.000)))/(1.00000+0.120000*exp((states[0]+12.0000)/10.0000))
algebraic[31] = 0.0500000*exp(-(states[0]+12.0000)/13.0000)
rates[12] = (4.00000*algebraic[16]*algebraic[2]+2.00000*algebraic[31]*states[13])-(algebraic[31]*states[12]+3.00000*algebraic[16]*states[12])
rates[13] = (3.00000*algebraic[16]*states[12]+3.00000*algebraic[31]*states[14])-(2.00000*algebraic[31]*states[13]+2.00000*algebraic[16]*states[13])
algebraic[13] = 1.00000-(states[32]+states[33]+states[31]+states[34])
algebraic[28] = 0.0223480*exp(0.0117600*states[0])
algebraic[33] = 0.0470020*exp(-0.0631000*states[0])
rates[32] = (algebraic[28]*algebraic[13]+constants[67]*states[33])-(algebraic[33]*states[32]+constants[68]*states[32])
algebraic[22] = (0.000152000*exp(-(states[0]-3.81000)/15.7500))/(0.00670830*exp(-(states[0]+132.050)/15.7500)+1.00000)
algebraic[32] = (0.000950000*exp((states[0]+132.050)/15.7500))/(0.0513350*exp((states[0]+132.050)/15.7500)+1.00000)
algebraic[37] = algebraic[22]/(algebraic[22]+algebraic[32])
algebraic[8] = 1.00000/((0.000152000*exp(-(states[0]+13.5000)/7.00000))/(0.0670830*exp(-(states[0]+33.5000)/7.00000)+1.00000)+(0.000950000*exp((states[0]+33.5000)/7.00000))/(0.0513350*exp((states[0]+33.5000)/7.00000)+1.00000))
rates[26] = (algebraic[37]-states[26])/algebraic[8]
algebraic[34] = 0.0908210*exp(0.0233910*(states[0]+5.00000))
algebraic[38] = 0.00649700*exp(-0.0326800*(states[0]+5.00000))
rates[31] = (algebraic[14]*states[33]+algebraic[38]*states[34])-(algebraic[29]*states[31]+algebraic[34]*states[31])
rates[34] = algebraic[34]*states[31]-algebraic[38]*states[34]
algebraic[36] = (constants[47]*states[2])/(constants[48]+states[2])
algebraic[40] = 13.0000*(1.00000-exp(-(power(states[0]+14.5000, 2.00000))/100.000))
rates[11] = (algebraic[16]*states[14]+constants[46]*states[15]+0.00100000*(algebraic[16]*states[16]-algebraic[40]*states[11]))-(4.00000*algebraic[31]*states[11]+algebraic[36]*states[11])
rates[14] = (2.00000*algebraic[16]*states[13]+4.00000*algebraic[31]*states[11]+0.0100000*(4.00000*constants[46]*algebraic[31]*states[15]-algebraic[16]*algebraic[36]*states[14])+0.00200000*(4.00000*algebraic[31]*states[16]-algebraic[40]*states[14])+4.00000*algebraic[31]*constants[46]*states[17])-(3.00000*algebraic[31]*states[14]+algebraic[16]*states[14]+1.00000*algebraic[36]*algebraic[40]*states[14])
rates[15] = (algebraic[36]*states[11]+0.00100000*(algebraic[16]*states[17]-algebraic[40]*states[15])+0.0100000*(algebraic[16]*algebraic[36]*states[14]-4.00000*algebraic[31]*algebraic[40]*states[15]))-constants[46]*states[15]
rates[16] = (0.00100000*(algebraic[40]*states[11]-algebraic[16]*states[16])+constants[46]*states[17]+0.00200000*(algebraic[40]*states[14]-4.00000*algebraic[31]*states[16]))-algebraic[36]*states[16]
rates[17] = (0.00100000*(algebraic[40]*states[15]-algebraic[16]*states[17])+algebraic[36]*states[16]+1.00000*algebraic[36]*algebraic[40]*states[14])-(4.00000*algebraic[31]*constants[46]*states[17]+constants[46]*states[17])
algebraic[35] = power(1.00000+(constants[20]*constants[22])/(power(constants[22]+states[3], 2.00000)), -1.00000)
algebraic[39] = constants[27]*(states[8]+states[9])*(states[3]-states[2])*states[5]
algebraic[42] = (states[2]-states[1])/constants[30]
rates[3] = algebraic[35]*(algebraic[42]-algebraic[39])
algebraic[43] = constants[45]*states[11]*(states[0]-constants[44])
rates[5] = -0.0400000*states[5]-((0.100000*algebraic[43])/constants[35])*exp(-(power(states[0]-5.00000, 2.00000))/648.000)
algebraic[44] = constants[29]*(states[4]-states[1])
algebraic[46] = (constants[31]*(power(states[1], 2.00000)))/(power(constants[32], 2.00000)+power(states[1], 2.00000))
algebraic[41] = (states[4]-states[3])/constants[28]
rates[4] = ((algebraic[46]-algebraic[44])*constants[1])/constants[3]-(algebraic[41]*constants[2])/constants[3]
algebraic[15] = power(1.00000+(constants[19]*constants[21])/(power(constants[21]+states[1], 2.00000)), -1.00000)
algebraic[48] = (constants[23]*states[1]*(constants[34]-states[7])+constants[25]*states[1]*(constants[33]-states[6]))-(constants[24]*states[7]+constants[26]*states[6])
rates[1] = algebraic[15]*((algebraic[44]+algebraic[42])-(algebraic[46]+algebraic[48]))
algebraic[45] = constants[50]*states[18]*states[19]*(states[0]-constants[49])
algebraic[47] = (constants[51]*(power(states[1], 2.00000)))/(power(constants[52], 2.00000)+power(states[1], 2.00000))
algebraic[49] = ((((((constants[53]*1.00000)/(power(constants[54], 3.00000)+power(constants[7], 3.00000)))*1.00000)/(constants[55]+constants[8]))*1.00000)/(1.00000+constants[56]*exp(((constants[57]-1.00000)*states[0]*constants[13])/(constants[11]*constants[12]))))*(exp((constants[57]*states[0]*constants[13])/(constants[11]*constants[12]))*(power(states[20], 3.00000))*constants[8]-exp(((constants[57]-1.00000)*states[0]*constants[13])/(constants[11]*constants[12]))*(power(constants[7], 3.00000))*states[1])
algebraic[50] = ((constants[11]*constants[12])/(2.00000*constants[13]))*log(constants[8]/states[1])
algebraic[51] = constants[58]*(states[0]-algebraic[50])
algebraic[30] = power(1.00000+(constants[19]*constants[21])/(power(constants[21]+states[2], 2.00000)), -1.00000)
rates[2] = algebraic[30]*((algebraic[39]*constants[2])/constants[4]-((algebraic[42]*constants[1])/constants[4]+((algebraic[43]+algebraic[51]+algebraic[47]+algebraic[45])-2.00000*algebraic[49])*((constants[5]*constants[0])/(2.00000*constants[4]*constants[13]))))
algebraic[52] = ((constants[11]*constants[12])/constants[13])*log(constants[7]/states[20])
algebraic[53] = constants[59]*(power(states[21], 3.00000))*states[22]*states[23]*(states[0]-algebraic[52])
algebraic[54] = constants[60]*(states[0]-algebraic[52])
algebraic[62] = 1.00000/(1.00000+0.124500*exp((-0.100000*states[0]*constants[13])/(constants[11]*constants[12]))+0.0365000*constants[75]*exp((-states[0]*constants[13])/(constants[11]*constants[12])))
algebraic[63] = (((constants[69]*algebraic[62]*1.00000)/(1.00000+power(constants[70]/states[20], 1.50000)))*constants[6])/(constants[6]+constants[71])
rates[20] = (-(algebraic[53]+algebraic[54]+3.00000*algebraic[63]+3.00000*algebraic[49])*constants[5]*constants[0])/(constants[1]*constants[13])
algebraic[55] = ((constants[11]*constants[12])/constants[13])*log(constants[6]/states[24])
algebraic[56] = constants[61]*(power(states[25], 6.50000))*states[26]*(states[0]-algebraic[55])
algebraic[59] = constants[64]*(constants[6]/(constants[6]+210.000))*((states[0]-algebraic[55])/(1.00000+exp(0.0896000*(states[0]-algebraic[55]))))
algebraic[60] = constants[65]*(power(states[30], 2.00000))*(states[0]-algebraic[55])
algebraic[57] = constants[62]*states[27]*states[28]*(states[0]-algebraic[55])
algebraic[58] = constants[63]*states[29]*(states[0]-algebraic[55])
algebraic[61] = constants[66]*states[31]*(states[0]-algebraic[55])
rates[24] = (-((algebraic[56]+algebraic[59]+algebraic[60]+algebraic[58]+algebraic[57]+algebraic[61])-2.00000*algebraic[63])*constants[5]*constants[0])/(constants[1]*constants[13])
algebraic[0] = custom_piecewise([greater_equal(voi , constants[14]) & less_equal(voi , constants[15]) & less_equal((voi-constants[14])-floor((voi-constants[14])/constants[16])*constants[16] , constants[17]), constants[18] , True, 0.00000])
algebraic[64] = power(states[2]/(constants[73]+states[2]), constants[74])
algebraic[65] = 1.00000*constants[72]*algebraic[64]*((states[0]*(power(constants[13], 2.00000)))/(constants[11]*constants[12]))*((constants[9]*exp((states[0]*constants[13])/(constants[11]*constants[12]))-constants[10])/(exp((states[0]*constants[13])/(constants[11]*constants[12]))-1.00000))
rates[0] = -(algebraic[43]+algebraic[45]+algebraic[47]+algebraic[49]+algebraic[51]+algebraic[53]+algebraic[54]+algebraic[63]+algebraic[56]+algebraic[59]+algebraic[60]+algebraic[57]+algebraic[58]+algebraic[61]+algebraic[65]+algebraic[0])
return(rates)
def computeAlgebraic(constants, states, voi):
algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
states = array(states)
voi = array(voi)
algebraic[1] = 1.00000-(states[10]+states[8]+states[9])
algebraic[3] = 1.00000/(1.00000+exp(-(states[0]+48.0000)/6.10000))
algebraic[17] = 0.100000+5.40000/(1.00000+exp((states[0]+100.000)/6.60000))
algebraic[4] = 1.00000/(1.00000+exp((states[0]+66.0000)/6.60000))
algebraic[18] = 8.00000+32.0000/(1.00000+exp((states[0]+65.0000)/5.00000))
algebraic[5] = (0.320000*(states[0]+47.1300))/(1.00000-exp(-0.100000*(states[0]+47.1300)))
algebraic[19] = 0.0800000*exp(-states[0]/11.0000)
algebraic[6] = custom_piecewise([less(states[0] , -40.0000), 0.135000*exp((80.0000+states[0])/-6.80000) , True, 0.00000])
algebraic[20] = custom_piecewise([less(states[0] , -40.0000), 3.56000*exp(0.0790000*states[0])+310000.*exp(0.350000*states[0]) , True, 1.00000/(0.130000*(1.00000+exp((states[0]+10.6600)/-11.1000)))])
algebraic[7] = custom_piecewise([less(states[0] , -40.0000), (-(127140.*exp(0.244400*states[0])+3.47400e-05*exp(-0.0439100*states[0]))*(states[0]+37.7800))/(1.00000+exp(0.311000*(states[0]+79.2300))) , True, 0.00000])
algebraic[21] = custom_piecewise([less(states[0] , -40.0000), (0.121200*exp(-0.0105200*states[0]))/(1.00000+exp(-0.137800*(states[0]+40.1400))) , True, (0.300000*exp(-2.53500e-07*states[0]))/(1.00000+exp(-0.100000*(states[0]+32.0000)))])
algebraic[9] = 0.180640*exp(0.0357700*(states[0]+30.0000))
algebraic[23] = 0.395600*exp(-0.0623700*(states[0]+30.0000))
algebraic[10] = 1.00000/(1.00000+exp(-(states[0]+22.5000)/7.70000))
algebraic[24] = 0.493000*exp(-0.0629000*states[0])+2.05800
algebraic[25] = 1200.00-170.000/(1.00000+exp((states[0]+45.2000)/5.70000))
algebraic[11] = 1.00000/(1.00000+exp((states[0]+45.2000)/5.70000))
algebraic[26] = 39.3000*exp(-0.0862000*states[0])+13.1700
algebraic[12] = (4.81333e-06*(states[0]+26.5000))/(1.00000-exp(-0.128000*(states[0]+26.5000)))
algebraic[27] = 9.53333e-05*exp(-0.0380000*(states[0]+26.5000))
algebraic[14] = 0.0137330*exp(0.0381980*states[0])
algebraic[29] = 6.89000e-05*exp(-0.0417800*states[0])
algebraic[2] = 1.00000-(states[11]+states[12]+states[12]+states[13]+states[14]+states[15]+states[16]+states[17])
algebraic[16] = (0.400000*exp((states[0]+12.0000)/10.0000)*((1.00000+0.700000*exp(-(power(states[0]+40.0000, 2.00000))/10.0000))-0.750000*exp(-(power(states[0]+20.0000, 2.00000))/400.000)))/(1.00000+0.120000*exp((states[0]+12.0000)/10.0000))
algebraic[31] = 0.0500000*exp(-(states[0]+12.0000)/13.0000)
algebraic[13] = 1.00000-(states[32]+states[33]+states[31]+states[34])
algebraic[28] = 0.0223480*exp(0.0117600*states[0])
algebraic[33] = 0.0470020*exp(-0.0631000*states[0])
algebraic[22] = (0.000152000*exp(-(states[0]-3.81000)/15.7500))/(0.00670830*exp(-(states[0]+132.050)/15.7500)+1.00000)
algebraic[32] = (0.000950000*exp((states[0]+132.050)/15.7500))/(0.0513350*exp((states[0]+132.050)/15.7500)+1.00000)
algebraic[37] = algebraic[22]/(algebraic[22]+algebraic[32])
algebraic[8] = 1.00000/((0.000152000*exp(-(states[0]+13.5000)/7.00000))/(0.0670830*exp(-(states[0]+33.5000)/7.00000)+1.00000)+(0.000950000*exp((states[0]+33.5000)/7.00000))/(0.0513350*exp((states[0]+33.5000)/7.00000)+1.00000))
algebraic[34] = 0.0908210*exp(0.0233910*(states[0]+5.00000))
algebraic[38] = 0.00649700*exp(-0.0326800*(states[0]+5.00000))
algebraic[36] = (constants[47]*states[2])/(constants[48]+states[2])
algebraic[40] = 13.0000*(1.00000-exp(-(power(states[0]+14.5000, 2.00000))/100.000))
algebraic[35] = power(1.00000+(constants[20]*constants[22])/(power(constants[22]+states[3], 2.00000)), -1.00000)
algebraic[39] = constants[27]*(states[8]+states[9])*(states[3]-states[2])*states[5]
algebraic[42] = (states[2]-states[1])/constants[30]
algebraic[43] = constants[45]*states[11]*(states[0]-constants[44])
algebraic[44] = constants[29]*(states[4]-states[1])
algebraic[46] = (constants[31]*(power(states[1], 2.00000)))/(power(constants[32], 2.00000)+power(states[1], 2.00000))
algebraic[41] = (states[4]-states[3])/constants[28]
algebraic[15] = power(1.00000+(constants[19]*constants[21])/(power(constants[21]+states[1], 2.00000)), -1.00000)
algebraic[48] = (constants[23]*states[1]*(constants[34]-states[7])+constants[25]*states[1]*(constants[33]-states[6]))-(constants[24]*states[7]+constants[26]*states[6])
algebraic[45] = constants[50]*states[18]*states[19]*(states[0]-constants[49])
algebraic[47] = (constants[51]*(power(states[1], 2.00000)))/(power(constants[52], 2.00000)+power(states[1], 2.00000))
algebraic[49] = ((((((constants[53]*1.00000)/(power(constants[54], 3.00000)+power(constants[7], 3.00000)))*1.00000)/(constants[55]+constants[8]))*1.00000)/(1.00000+constants[56]*exp(((constants[57]-1.00000)*states[0]*constants[13])/(constants[11]*constants[12]))))*(exp((constants[57]*states[0]*constants[13])/(constants[11]*constants[12]))*(power(states[20], 3.00000))*constants[8]-exp(((constants[57]-1.00000)*states[0]*constants[13])/(constants[11]*constants[12]))*(power(constants[7], 3.00000))*states[1])
algebraic[50] = ((constants[11]*constants[12])/(2.00000*constants[13]))*log(constants[8]/states[1])
algebraic[51] = constants[58]*(states[0]-algebraic[50])
algebraic[30] = power(1.00000+(constants[19]*constants[21])/(power(constants[21]+states[2], 2.00000)), -1.00000)
algebraic[52] = ((constants[11]*constants[12])/constants[13])*log(constants[7]/states[20])
algebraic[53] = constants[59]*(power(states[21], 3.00000))*states[22]*states[23]*(states[0]-algebraic[52])
algebraic[54] = constants[60]*(states[0]-algebraic[52])
algebraic[62] = 1.00000/(1.00000+0.124500*exp((-0.100000*states[0]*constants[13])/(constants[11]*constants[12]))+0.0365000*constants[75]*exp((-states[0]*constants[13])/(constants[11]*constants[12])))
algebraic[63] = (((constants[69]*algebraic[62]*1.00000)/(1.00000+power(constants[70]/states[20], 1.50000)))*constants[6])/(constants[6]+constants[71])
algebraic[55] = ((constants[11]*constants[12])/constants[13])*log(constants[6]/states[24])
algebraic[56] = constants[61]*(power(states[25], 6.50000))*states[26]*(states[0]-algebraic[55])
algebraic[59] = constants[64]*(constants[6]/(constants[6]+210.000))*((states[0]-algebraic[55])/(1.00000+exp(0.0896000*(states[0]-algebraic[55]))))
algebraic[60] = constants[65]*(power(states[30], 2.00000))*(states[0]-algebraic[55])
algebraic[57] = constants[62]*states[27]*states[28]*(states[0]-algebraic[55])
algebraic[58] = constants[63]*states[29]*(states[0]-algebraic[55])
algebraic[61] = constants[66]*states[31]*(states[0]-algebraic[55])
algebraic[0] = custom_piecewise([greater_equal(voi , constants[14]) & less_equal(voi , constants[15]) & less_equal((voi-constants[14])-floor((voi-constants[14])/constants[16])*constants[16] , constants[17]), constants[18] , True, 0.00000])
algebraic[64] = power(states[2]/(constants[73]+states[2]), constants[74])
algebraic[65] = 1.00000*constants[72]*algebraic[64]*((states[0]*(power(constants[13], 2.00000)))/(constants[11]*constants[12]))*((constants[9]*exp((states[0]*constants[13])/(constants[11]*constants[12]))-constants[10])/(exp((states[0]*constants[13])/(constants[11]*constants[12]))-1.00000))
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)