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 = 71
sizeStates = 41
sizeConstants = 73
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] = "R in component membrane (joule_per_mole_kelvin)"
legend_constants[10] = "T in component membrane (kelvin)"
legend_constants[11] = "F in component membrane (coulomb_per_millimole)"
legend_algebraic[0] = "i_stim in component membrane (picoA_per_picoF)"
legend_algebraic[46] = "i_CaL in component L_type_calcium_current (picoA_per_picoF)"
legend_algebraic[48] = "i_pCa in component calcium_pump_current (picoA_per_picoF)"
legend_algebraic[50] = "i_NaCa in component sodium_calcium_exchange_current (picoA_per_picoF)"
legend_algebraic[54] = "i_Cab in component calcium_background_current (picoA_per_picoF)"
legend_algebraic[57] = "i_Na in component fast_sodium_current (picoA_per_picoF)"
legend_algebraic[58] = "i_Nab in component sodium_background_current (picoA_per_picoF)"
legend_algebraic[68] = "i_NaK in component sodium_potassium_pump_current (picoA_per_picoF)"
legend_algebraic[60] = "i_Kto_f in component fast_transient_outward_potassium_current (picoA_per_picoF)"
legend_algebraic[61] = "i_Kto_s in component slow_transient_outward_potassium_current (picoA_per_picoF)"
legend_algebraic[62] = "i_K1 in component time_independent_potassium_current (picoA_per_picoF)"
legend_algebraic[63] = "i_Ks in component slow_delayed_rectifier_potassium_current (picoA_per_picoF)"
legend_algebraic[64] = "i_Kur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (picoA_per_picoF)"
legend_algebraic[65] = "i_Kss in component non_inactivating_steady_state_potassium_current (picoA_per_picoF)"
legend_algebraic[70] = "i_ClCa in component calcium_activated_chloride_current (picoA_per_picoF)"
legend_algebraic[66] = "i_Kr in component rapid_delayed_rectifier_potassium_current (picoA_per_picoF)"
legend_constants[12] = "stim_start in component membrane (millisecond)"
legend_constants[13] = "stim_end in component membrane (millisecond)"
legend_constants[14] = "stim_period in component membrane (millisecond)"
legend_constants[15] = "stim_duration in component membrane (millisecond)"
legend_constants[16] = "stim_amplitude in component membrane (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[11] = "Bi in component calcium_concentration (dimensionless)"
legend_algebraic[24] = "Bss in component calcium_concentration (dimensionless)"
legend_algebraic[29] = "BJSR in component calcium_concentration (dimensionless)"
legend_constants[17] = "CMDN_tot in component calcium_concentration (micromolar)"
legend_constants[18] = "CSQN_tot in component calcium_concentration (micromolar)"
legend_constants[19] = "Km_CMDN in component calcium_concentration (micromolar)"
legend_constants[20] = "Km_CSQN in component calcium_concentration (micromolar)"
legend_algebraic[40] = "J_leak in component calcium_fluxes (micromolar_per_millisecond)"
legend_algebraic[33] = "J_rel in component calcium_fluxes (micromolar_per_millisecond)"
legend_algebraic[42] = "J_up in component calcium_fluxes (micromolar_per_millisecond)"
legend_algebraic[36] = "J_tr in component calcium_fluxes (micromolar_per_millisecond)"
legend_algebraic[44] = "J_trpn in component calcium_fluxes (micromolar_per_millisecond)"
legend_algebraic[38] = "J_xfer in component calcium_fluxes (micromolar_per_millisecond)"
legend_constants[21] = "k_plus_htrpn in component calcium_fluxes (per_micromolar_millisecond)"
legend_constants[22] = "k_minus_htrpn in component calcium_fluxes (per_millisecond)"
legend_constants[23] = "k_plus_ltrpn in component calcium_fluxes (per_micromolar_millisecond)"
legend_constants[24] = "k_minus_ltrpn in component calcium_fluxes (per_millisecond)"
legend_states[5] = "P_RyR in component calcium_fluxes (dimensionless)"
legend_constants[25] = "v1 in component calcium_fluxes (per_millisecond)"
legend_constants[26] = "tau_tr in component calcium_fluxes (millisecond)"
legend_constants[27] = "v2 in component calcium_fluxes (per_millisecond)"
legend_constants[28] = "tau_xfer in component calcium_fluxes (millisecond)"
legend_constants[29] = "v3 in component calcium_fluxes (micromolar_per_millisecond)"
legend_constants[30] = "Km_up in component calcium_fluxes (micromolar)"
legend_constants[31] = "LTRPN_tot in component calcium_buffering (micromolar)"
legend_constants[32] = "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[33] = "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[34] = "k_plus_a in component ryanodine_receptors (micromolar4_per_millisecond)"
legend_constants[35] = "k_minus_a in component ryanodine_receptors (per_millisecond)"
legend_constants[36] = "k_plus_b in component ryanodine_receptors (micromolar3_per_millisecond)"
legend_constants[37] = "k_minus_b in component ryanodine_receptors (per_millisecond)"
legend_constants[38] = "k_plus_c in component ryanodine_receptors (per_millisecond)"
legend_constants[39] = "k_minus_c in component ryanodine_receptors (per_millisecond)"
legend_constants[40] = "m in component ryanodine_receptors (dimensionless)"
legend_constants[41] = "n in component ryanodine_receptors (dimensionless)"
legend_constants[42] = "E_CaL in component L_type_calcium_current (millivolt)"
legend_constants[43] = "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[12] = "alpha in component L_type_calcium_current (per_millisecond)"
legend_algebraic[25] = "beta in component L_type_calcium_current (per_millisecond)"
legend_algebraic[30] = "gamma in component L_type_calcium_current (per_millisecond)"
legend_algebraic[34] = "Kpcf in component L_type_calcium_current (per_millisecond)"
legend_constants[44] = "Kpcb in component L_type_calcium_current (per_millisecond)"
legend_constants[45] = "Kpc_max in component L_type_calcium_current (per_millisecond)"
legend_constants[46] = "Kpc_half in component L_type_calcium_current (micromolar)"
legend_constants[47] = "i_pCa_max in component calcium_pump_current (picoA_per_picoF)"
legend_constants[48] = "Km_pCa in component calcium_pump_current (micromolar)"
legend_constants[49] = "k_NaCa in component sodium_calcium_exchange_current (picoA_per_picoF)"
legend_constants[50] = "K_mNa in component sodium_calcium_exchange_current (micromolar)"
legend_constants[51] = "K_mCa in component sodium_calcium_exchange_current (micromolar)"
legend_constants[52] = "k_sat in component sodium_calcium_exchange_current (dimensionless)"
legend_constants[53] = "eta in component sodium_calcium_exchange_current (dimensionless)"
legend_states[18] = "Nai in component sodium_concentration (micromolar)"
legend_constants[54] = "g_Cab in component calcium_background_current (milliS_per_microF)"
legend_algebraic[52] = "E_CaN in component calcium_background_current (millivolt)"
legend_algebraic[56] = "E_Na in component fast_sodium_current (millivolt)"
legend_constants[55] = "g_Na in component fast_sodium_current (milliS_per_microF)"
legend_states[19] = "O_Na in component fast_sodium_current (dimensionless)"
legend_states[20] = "C_Na1 in component fast_sodium_current (dimensionless)"
legend_states[21] = "C_Na2 in component fast_sodium_current (dimensionless)"
legend_algebraic[3] = "C_Na3 in component fast_sodium_current (dimensionless)"
legend_states[22] = "I1_Na in component fast_sodium_current (dimensionless)"
legend_states[23] = "I2_Na in component fast_sodium_current (dimensionless)"
legend_states[24] = "IF_Na in component fast_sodium_current (dimensionless)"
legend_states[25] = "IC_Na2 in component fast_sodium_current (dimensionless)"
legend_states[26] = "IC_Na3 in component fast_sodium_current (dimensionless)"
legend_algebraic[13] = "alpha_Na11 in component fast_sodium_current (per_millisecond)"
legend_algebraic[35] = "beta_Na11 in component fast_sodium_current (per_millisecond)"
legend_algebraic[26] = "alpha_Na12 in component fast_sodium_current (per_millisecond)"
legend_algebraic[37] = "beta_Na12 in component fast_sodium_current (per_millisecond)"
legend_algebraic[31] = "alpha_Na13 in component fast_sodium_current (per_millisecond)"
legend_algebraic[39] = "beta_Na13 in component fast_sodium_current (per_millisecond)"
legend_algebraic[41] = "alpha_Na3 in component fast_sodium_current (per_millisecond)"
legend_algebraic[43] = "beta_Na3 in component fast_sodium_current (per_millisecond)"
legend_algebraic[45] = "alpha_Na2 in component fast_sodium_current (per_millisecond)"
legend_algebraic[47] = "beta_Na2 in component fast_sodium_current (per_millisecond)"
legend_algebraic[49] = "alpha_Na4 in component fast_sodium_current (per_millisecond)"
legend_algebraic[51] = "beta_Na4 in component fast_sodium_current (per_millisecond)"
legend_algebraic[53] = "alpha_Na5 in component fast_sodium_current (per_millisecond)"
legend_algebraic[55] = "beta_Na5 in component fast_sodium_current (per_millisecond)"
legend_states[27] = "Ki in component potassium_concentration (micromolar)"
legend_constants[56] = "g_Nab in component sodium_background_current (milliS_per_microF)"
legend_algebraic[59] = "E_K in component fast_transient_outward_potassium_current (millivolt)"
legend_constants[57] = "g_Kto_f in component fast_transient_outward_potassium_current (milliS_per_microF)"
legend_states[28] = "ato_f in component fast_transient_outward_potassium_current (dimensionless)"
legend_states[29] = "ito_f in component fast_transient_outward_potassium_current (dimensionless)"
legend_algebraic[4] = "alpha_a in component fast_transient_outward_potassium_current (per_millisecond)"
legend_algebraic[14] = "beta_a in component fast_transient_outward_potassium_current (per_millisecond)"
legend_algebraic[5] = "alpha_i in component fast_transient_outward_potassium_current (per_millisecond)"
legend_algebraic[15] = "beta_i in component fast_transient_outward_potassium_current (per_millisecond)"
legend_algebraic[6] = "ass in component slow_transient_outward_potassium_current (dimensionless)"
legend_algebraic[7] = "iss in component slow_transient_outward_potassium_current (dimensionless)"
legend_constants[58] = "g_Kto_s in component slow_transient_outward_potassium_current (milliS_per_microF)"
legend_states[30] = "ato_s in component slow_transient_outward_potassium_current (dimensionless)"
legend_states[31] = "ito_s in component slow_transient_outward_potassium_current (dimensionless)"
legend_algebraic[16] = "tau_ta_s in component slow_transient_outward_potassium_current (millisecond)"
legend_algebraic[17] = "tau_ti_s in component slow_transient_outward_potassium_current (millisecond)"
legend_constants[59] = "g_Ks in component slow_delayed_rectifier_potassium_current (milliS_per_microF)"
legend_states[32] = "nKs in component slow_delayed_rectifier_potassium_current (dimensionless)"
legend_algebraic[8] = "alpha_n in component slow_delayed_rectifier_potassium_current (per_millisecond)"
legend_algebraic[18] = "beta_n in component slow_delayed_rectifier_potassium_current (per_millisecond)"
legend_constants[60] = "g_Kur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (milliS_per_microF)"
legend_states[33] = "aur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)"
legend_states[34] = "iur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)"
legend_algebraic[19] = "tau_aur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (millisecond)"
legend_algebraic[20] = "tau_iur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (millisecond)"
legend_constants[61] = "g_Kss in component non_inactivating_steady_state_potassium_current (milliS_per_microF)"
legend_states[35] = "aKss in component non_inactivating_steady_state_potassium_current (dimensionless)"
legend_states[36] = "iKss in component non_inactivating_steady_state_potassium_current (dimensionless)"
legend_algebraic[21] = "tau_Kss in component non_inactivating_steady_state_potassium_current (millisecond)"
legend_constants[62] = "g_Kr in component rapid_delayed_rectifier_potassium_current (milliS_per_microF)"
legend_states[37] = "O_K in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_states[38] = "C_K1 in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_states[39] = "C_K2 in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_algebraic[9] = "C_K0 in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_states[40] = "I_K in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_algebraic[22] = "alpha_a0 in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_algebraic[27] = "beta_a0 in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_constants[63] = "kb in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_constants[64] = "kf in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_algebraic[10] = "alpha_a1 in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_algebraic[23] = "beta_a1 in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_algebraic[28] = "alpha_i in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_algebraic[32] = "beta_i in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
legend_constants[65] = "i_NaK_max in component sodium_potassium_pump_current (picoA_per_picoF)"
legend_constants[66] = "Km_Nai in component sodium_potassium_pump_current (micromolar)"
legend_constants[67] = "Km_Ko in component sodium_potassium_pump_current (micromolar)"
legend_algebraic[67] = "f_NaK in component sodium_potassium_pump_current (dimensionless)"
legend_constants[71] = "sigma in component sodium_potassium_pump_current (dimensionless)"
legend_constants[68] = "g_ClCa in component calcium_activated_chloride_current (milliS_per_microF)"
legend_algebraic[69] = "O_ClCa in component calcium_activated_chloride_current (dimensionless)"
legend_constants[69] = "E_Cl in component calcium_activated_chloride_current (millivolt)"
legend_constants[70] = "Km_Cl in component calcium_activated_chloride_current (micromolar)"
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 Nai in component sodium_concentration (micromolar)"
legend_rates[21] = "d/dt C_Na2 in component fast_sodium_current (dimensionless)"
legend_rates[20] = "d/dt C_Na1 in component fast_sodium_current (dimensionless)"
legend_rates[19] = "d/dt O_Na in component fast_sodium_current (dimensionless)"
legend_rates[24] = "d/dt IF_Na in component fast_sodium_current (dimensionless)"
legend_rates[22] = "d/dt I1_Na in component fast_sodium_current (dimensionless)"
legend_rates[23] = "d/dt I2_Na in component fast_sodium_current (dimensionless)"
legend_rates[25] = "d/dt IC_Na2 in component fast_sodium_current (dimensionless)"
legend_rates[26] = "d/dt IC_Na3 in component fast_sodium_current (dimensionless)"
legend_rates[27] = "d/dt Ki in component potassium_concentration (micromolar)"
legend_rates[28] = "d/dt ato_f in component fast_transient_outward_potassium_current (dimensionless)"
legend_rates[29] = "d/dt ito_f in component fast_transient_outward_potassium_current (dimensionless)"
legend_rates[30] = "d/dt ato_s in component slow_transient_outward_potassium_current (dimensionless)"
legend_rates[31] = "d/dt ito_s in component slow_transient_outward_potassium_current (dimensionless)"
legend_rates[32] = "d/dt nKs in component slow_delayed_rectifier_potassium_current (dimensionless)"
legend_rates[33] = "d/dt aur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)"
legend_rates[34] = "d/dt iur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)"
legend_rates[35] = "d/dt aKss in component non_inactivating_steady_state_potassium_current (dimensionless)"
legend_rates[36] = "d/dt iKss in component non_inactivating_steady_state_potassium_current (dimensionless)"
legend_rates[39] = "d/dt C_K2 in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_rates[38] = "d/dt C_K1 in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_rates[37] = "d/dt O_K in component rapid_delayed_rectifier_potassium_current (dimensionless)"
legend_rates[40] = "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] = -82.4202
constants[0] = 1
constants[1] = 25.84e-6
constants[2] = 0.12e-6
constants[3] = 2.098e-6
constants[4] = 1.485e-9
constants[5] = 1.534e-4
constants[6] = 5400
constants[7] = 140000
constants[8] = 1800
constants[9] = 8.314
constants[10] = 298
constants[11] = 96.5
constants[12] = 20
constants[13] = 100000
constants[14] = 71.43
constants[15] = 0.5
constants[16] = -80
states[1] = 0.115001
states[2] = 0.115001
states[3] = 1299.5
states[4] = 1299.5
constants[17] = 50
constants[18] = 15000
constants[19] = 0.238
constants[20] = 800
constants[21] = 0.00237
constants[22] = 3.2e-5
constants[23] = 0.0327
constants[24] = 0.0196
states[5] = 0
constants[25] = 4.5
constants[26] = 20
constants[27] = 1.74e-5
constants[28] = 8
constants[29] = 0.45
constants[30] = 0.5
constants[31] = 70
constants[32] = 140
states[6] = 11.2684
states[7] = 125.29
constants[33] = 7
states[8] = 0.149102e-4
states[9] = 0.951726e-10
states[10] = 0.16774e-3
constants[34] = 0.006075
constants[35] = 0.07125
constants[36] = 0.00405
constants[37] = 0.965
constants[38] = 0.009
constants[39] = 0.0008
constants[40] = 3
constants[41] = 4
constants[42] = 63
constants[43] = 0.1729
states[11] = 0.930308e-18
states[12] = 0.124216e-3
states[13] = 0.578679e-8
states[14] = 0.119816e-12
states[15] = 0.497923e-18
states[16] = 0.345847e-13
states[17] = 0.185106e-13
constants[44] = 0.0005
constants[45] = 0.23324
constants[46] = 20
constants[47] = 1
constants[48] = 0.5
constants[49] = 292.8
constants[50] = 87500
constants[51] = 1380
constants[52] = 0.1
constants[53] = 0.35
states[18] = 14237.1
constants[54] = 0.000367
constants[55] = 13
states[19] = 0.713483e-6
states[20] = 0.279132e-3
states[21] = 0.020752
states[22] = 0.673345e-6
states[23] = 0.155787e-8
states[24] = 0.153176e-3
states[25] = 0.0113879
states[26] = 0.34278
states[27] = 143720
constants[56] = 0.0026
constants[57] = 0.0798
states[28] = 0.265563e-2
states[29] = 0.999977
constants[58] = 0.0629
states[30] = 0.417069e-3
states[31] = 0.998543
constants[59] = 0.00575
states[32] = 0.262753e-3
constants[60] = 0.0975
states[33] = 0.417069e-3
states[34] = 0.998543
constants[61] = 0.0324
states[35] = 0.417069e-3
states[36] = 1
constants[62] = 0.078
states[37] = 0.175298e-3
states[38] = 0.992513e-3
states[39] = 0.641229e-3
states[40] = 0.319129e-4
constants[63] = 0.036778
constants[64] = 0.023761
constants[65] = 0.88
constants[66] = 21000
constants[67] = 1500
constants[68] = 10
constants[69] = -40
constants[70] = 10
constants[71] = (1.00000/7.00000)*(exp(constants[7]/67300.0)-1.00000)
constants[72] = 0.00000
return (states, constants)
def computeRates(voi, states, constants):
rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
rates[36] = constants[72]
rates[6] = constants[23]*states[1]*(constants[31]-states[6])-constants[24]*states[6]
rates[7] = constants[21]*states[1]*(constants[32]-states[7])-constants[22]*states[7]
rates[9] = constants[36]*(power(states[2], constants[40]))*states[8]-constants[37]*states[9]
rates[10] = constants[38]*states[8]-constants[39]*states[10]
algebraic[1] = 1.00000-(states[10]+states[8]+states[9])
rates[8] = (constants[34]*(power(states[2], constants[41]))*algebraic[1]+constants[37]*states[9]+constants[39]*states[10])-(constants[35]*states[8]+constants[36]*(power(states[2], constants[40]))*states[8]+constants[38]*states[8])
algebraic[4] = 0.180640*exp(0.0357700*(states[0]+30.0000))
algebraic[14] = 0.395600*exp(-0.0623700*(states[0]+30.0000))
rates[28] = algebraic[4]*(1.00000-states[28])-algebraic[14]*states[28]
algebraic[5] = (0.000152000*exp(-(states[0]+13.5000)/7.00000))/(0.00670830*exp(-(states[0]+33.5000)/7.00000)+1.00000)
algebraic[15] = (0.000950000*exp((states[0]+33.5000)/7.00000))/(0.0513350*exp((states[0]+33.5000)/7.00000)+1.00000)
rates[29] = algebraic[5]*(1.00000-states[29])-algebraic[15]*states[29]
algebraic[6] = 1.00000/(1.00000+exp(-(states[0]+22.5000)/7.70000))
algebraic[16] = 0.493000*exp(-0.0629000*states[0])+2.05800
rates[30] = (algebraic[6]-states[30])/algebraic[16]
algebraic[7] = 1.00000/(1.00000+exp((states[0]+45.2000)/5.70000))
algebraic[17] = 270.000+1050.00/(1.00000+exp((states[0]+45.2000)/5.70000))
rates[31] = (algebraic[7]-states[31])/algebraic[17]
algebraic[8] = (4.81333e-06*(states[0]+26.5000))/(1.00000-exp(-0.128000*(states[0]+26.5000)))
algebraic[18] = 9.53333e-05*exp(-0.0380000*(states[0]+26.5000))
rates[32] = algebraic[8]*(1.00000-states[32])-algebraic[18]*states[32]
algebraic[19] = 0.493000*exp(-0.0629000*states[0])+2.05800
rates[33] = (algebraic[6]-states[33])/algebraic[19]
algebraic[20] = 1200.00-170.000/(1.00000+exp((states[0]+45.2000)/5.70000))
rates[34] = (algebraic[7]-states[34])/algebraic[20]
algebraic[21] = 39.3000*exp(-0.0862000*states[0])+13.1700
rates[35] = (algebraic[6]-states[35])/algebraic[21]
algebraic[10] = 0.0137330*exp(0.0381980*states[0])
algebraic[23] = 6.89000e-05*exp(-0.0417800*states[0])
rates[39] = (constants[64]*states[38]+algebraic[23]*states[37])-(constants[63]*states[39]+algebraic[10]*states[39])
algebraic[2] = 1.00000-(states[11]+states[12]+states[13]+states[14]+states[15]+states[16]+states[17])
algebraic[12] = (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[25] = 0.0500000*exp(-(states[0]+12.0000)/13.0000)
rates[12] = (4.00000*algebraic[12]*algebraic[2]+2.00000*algebraic[25]*states[13])-(algebraic[25]*states[12]+3.00000*algebraic[12]*states[12])
rates[13] = (3.00000*algebraic[12]*states[12]+3.00000*algebraic[25]*states[14])-(2.00000*algebraic[25]*states[13]+2.00000*algebraic[12]*states[13])
algebraic[9] = 1.00000-(states[38]+states[39]+states[37]+states[40])
algebraic[22] = 0.0223480*exp(0.0117600*states[0])
algebraic[27] = 0.0470020*exp(-0.0631000*states[0])
rates[38] = (algebraic[22]*algebraic[9]+constants[63]*states[39])-(algebraic[27]*states[38]+constants[64]*states[38])
algebraic[28] = 0.0908210*exp(0.0233910*(states[0]+5.00000))
algebraic[32] = 0.00649700*exp(-0.0326800*(states[0]+5.00000))
rates[37] = (algebraic[10]*states[39]+algebraic[32]*states[40])-(algebraic[23]*states[37]+algebraic[28]*states[37])
rates[40] = algebraic[28]*states[37]-algebraic[32]*states[40]
algebraic[30] = (constants[45]*states[2])/(constants[46]+states[2])
algebraic[34] = 13.0000*(1.00000-exp(-(power(states[0]+14.5000, 2.00000))/100.000))
rates[11] = (algebraic[12]*states[14]+constants[44]*states[15]+0.00100000*(algebraic[12]*states[16]-algebraic[34]*states[11]))-(4.00000*algebraic[25]*states[11]+algebraic[30]*states[11])
rates[14] = (2.00000*algebraic[12]*states[13]+4.00000*algebraic[25]*states[11]+0.0100000*(4.00000*constants[44]*algebraic[25]*states[15]-algebraic[12]*algebraic[30]*states[14])+0.00200000*(4.00000*algebraic[25]*states[16]-algebraic[34]*states[14])+4.00000*algebraic[25]*constants[44]*states[17])-(3.00000*algebraic[25]*states[14]+algebraic[12]*states[14]+1.00000*algebraic[30]*algebraic[34]*states[14])
rates[15] = (algebraic[30]*states[11]+0.00100000*(algebraic[12]*states[17]-algebraic[34]*states[15])+0.0100000*(algebraic[12]*algebraic[30]*states[14]-4.00000*algebraic[25]*algebraic[34]*states[15]))-constants[44]*states[15]
rates[16] = (0.00100000*(algebraic[34]*states[11]-algebraic[12]*states[16])+constants[44]*states[17]+0.00200000*(algebraic[34]*states[14]-4.00000*algebraic[25]*states[16]))-algebraic[30]*states[16]
rates[17] = (0.00100000*(algebraic[34]*states[15]-algebraic[12]*states[17])+algebraic[30]*states[16]+1.00000*algebraic[30]*algebraic[34]*states[14])-(4.00000*algebraic[25]*constants[44]*states[17]+constants[44]*states[17])
algebraic[29] = power(1.00000+(constants[18]*constants[20])/(power(constants[20]+states[3], 2.00000)), -1.00000)
algebraic[33] = constants[25]*(states[8]+states[9])*(states[3]-states[2])*states[5]
algebraic[36] = (states[4]-states[3])/constants[26]
rates[3] = algebraic[29]*(algebraic[36]-algebraic[33])
algebraic[40] = constants[27]*(states[4]-states[1])
algebraic[42] = (constants[29]*(power(states[1], 2.00000)))/(power(constants[30], 2.00000)+power(states[1], 2.00000))
rates[4] = ((algebraic[42]-algebraic[40])*constants[1])/constants[3]-(algebraic[36]*constants[2])/constants[3]
algebraic[3] = 1.00000-(states[19]+states[20]+states[21]+states[24]+states[22]+states[23]+states[25]+states[26])
algebraic[13] = 3.80200/(0.102700*exp(-(states[0]+2.50000)/17.0000)+0.200000*exp(-(states[0]+2.50000)/150.000))
algebraic[35] = 0.191700*exp(-(states[0]+2.50000)/20.3000)
algebraic[26] = 3.80200/(0.102700*exp(-(states[0]+2.50000)/15.0000)+0.230000*exp(-(states[0]+2.50000)/150.000))
algebraic[37] = 0.200000*exp(-(states[0]-2.50000)/20.3000)
algebraic[41] = 7.00000e-07*exp(-(states[0]+7.00000)/7.70000)
algebraic[43] = 0.00840000+2.00000e-05*(states[0]+7.00000)
rates[21] = (algebraic[13]*algebraic[3]+algebraic[37]*states[20]+algebraic[41]*states[25])-(algebraic[35]*states[21]+algebraic[26]*states[21]+algebraic[43]*states[21])
algebraic[31] = 3.80200/(0.102700*exp(-(states[0]+2.50000)/12.0000)+0.250000*exp(-(states[0]+2.50000)/150.000))
algebraic[39] = 0.220000*exp(-(states[0]-7.50000)/20.3000)
rates[20] = (algebraic[26]*states[21]+algebraic[39]*states[19]+algebraic[41]*states[24])-(algebraic[37]*states[20]+algebraic[31]*states[20]+algebraic[43]*states[20])
rates[25] = (algebraic[13]*states[26]+algebraic[37]*states[24]+algebraic[43]*states[21])-(algebraic[35]*states[25]+algebraic[26]*states[25]+algebraic[41]*states[25])
rates[26] = (algebraic[35]*states[25]+algebraic[43]*algebraic[3])-(algebraic[13]*states[26]+algebraic[41]*states[26])
algebraic[46] = constants[43]*states[11]*(states[0]-constants[42])
algebraic[24] = power(1.00000+(constants[17]*constants[19])/(power(constants[19]+states[2], 2.00000)), -1.00000)
algebraic[38] = (states[2]-states[1])/constants[28]
rates[2] = algebraic[24]*((algebraic[33]*constants[2])/constants[4]-((algebraic[38]*constants[1])/constants[4]+(algebraic[46]*constants[5]*constants[0])/(2.00000*constants[4]*constants[11])))
rates[5] = -0.0400000*states[5]-((0.100000*algebraic[46])/constants[33])*exp(-(power(states[0]-5.00000, 2.00000))/648.000)
algebraic[45] = 1.00000/(0.188495*exp(-(states[0]+7.00000)/16.6000)+0.393956)
algebraic[47] = (algebraic[31]*algebraic[45]*algebraic[41])/(algebraic[39]*algebraic[43])
rates[19] = (algebraic[31]*states[20]+algebraic[47]*states[24])-(algebraic[39]*states[19]+algebraic[45]*states[19])
algebraic[49] = algebraic[45]/1000.00
algebraic[51] = algebraic[41]
rates[24] = (algebraic[45]*states[19]+algebraic[43]*states[20]+algebraic[51]*states[22]+algebraic[26]*states[25])-(algebraic[47]*states[24]+algebraic[41]*states[24]+algebraic[49]*states[24]+algebraic[37]*states[24])
algebraic[48] = (constants[47]*(power(states[1], 2.00000)))/(power(constants[48], 2.00000)+power(states[1], 2.00000))
algebraic[50] = ((((((constants[49]*1.00000)/(power(constants[50], 3.00000)+power(constants[7], 3.00000)))*1.00000)/(constants[51]+constants[8]))*1.00000)/(1.00000+constants[52]*exp(((constants[53]-1.00000)*states[0]*constants[11])/(constants[9]*constants[10]))))*(exp((constants[53]*states[0]*constants[11])/(constants[9]*constants[10]))*(power(states[18], 3.00000))*constants[8]-exp(((constants[53]-1.00000)*states[0]*constants[11])/(constants[9]*constants[10]))*(power(constants[7], 3.00000))*states[1])
algebraic[52] = ((constants[9]*constants[10])/(2.00000*constants[11]))*log(constants[8]/states[1])
algebraic[54] = constants[54]*(states[0]-algebraic[52])
algebraic[11] = power(1.00000+(constants[17]*constants[19])/(power(constants[19]+states[1], 2.00000)), -1.00000)
algebraic[44] = (constants[21]*states[1]*(constants[32]-states[7])+constants[23]*states[1]*(constants[31]-states[6]))-(constants[22]*states[7]+constants[24]*states[6])
rates[1] = algebraic[11]*((algebraic[40]+algebraic[38])-(algebraic[42]+algebraic[44]+(((algebraic[54]+algebraic[48])-2.00000*algebraic[50])*constants[5]*constants[0])/(2.00000*constants[1]*constants[11])))
algebraic[53] = algebraic[45]/95000.0
algebraic[55] = algebraic[41]/50.0000
rates[22] = (algebraic[49]*states[24]+algebraic[55]*states[23])-(algebraic[51]*states[22]+algebraic[53]*states[22])
rates[23] = algebraic[53]*states[22]-algebraic[55]*states[23]
algebraic[56] = ((constants[9]*constants[10])/constants[11])*log((0.900000*constants[7]+0.100000*constants[6])/(0.900000*states[18]+0.100000*states[27]))
algebraic[57] = constants[55]*states[19]*(states[0]-algebraic[56])
algebraic[58] = constants[56]*(states[0]-algebraic[56])
algebraic[67] = 1.00000/(1.00000+0.124500*exp((-0.100000*states[0]*constants[11])/(constants[9]*constants[10]))+0.0365000*constants[71]*exp((-states[0]*constants[11])/(constants[9]*constants[10])))
algebraic[68] = (((constants[65]*algebraic[67]*1.00000)/(1.00000+power(constants[66]/states[18], 1.50000)))*constants[6])/(constants[6]+constants[67])
rates[18] = (-(algebraic[57]+algebraic[58]+3.00000*algebraic[68]+3.00000*algebraic[50])*constants[5]*constants[0])/(constants[1]*constants[11])
algebraic[59] = ((constants[9]*constants[10])/constants[11])*log(constants[6]/states[27])
algebraic[60] = constants[57]*(power(states[28], 3.00000))*states[29]*(states[0]-algebraic[59])
algebraic[61] = constants[58]*states[30]*states[31]*(states[0]-algebraic[59])
algebraic[62] = (((0.293800*constants[6])/(constants[6]+210.000))*(states[0]-algebraic[59]))/(1.00000+exp(0.0896000*(states[0]-algebraic[59])))
algebraic[63] = constants[59]*(power(states[32], 2.00000))*(states[0]-algebraic[59])
algebraic[64] = constants[60]*states[33]*states[34]*(states[0]-algebraic[59])
algebraic[65] = constants[61]*states[35]*states[36]*(states[0]-algebraic[59])
algebraic[66] = constants[62]*states[37]*(states[0]-((constants[9]*constants[10])/constants[11])*log((0.980000*constants[6]+0.0200000*constants[7])/(0.980000*states[27]+0.0200000*states[18])))
rates[27] = (-((algebraic[60]+algebraic[61]+algebraic[62]+algebraic[63]+algebraic[65]+algebraic[64]+algebraic[66])-2.00000*algebraic[68])*constants[5]*constants[0])/(constants[1]*constants[11])
algebraic[0] = custom_piecewise([greater_equal(voi , constants[12]) & less_equal(voi , constants[13]) & less_equal((voi-constants[12])-floor((voi-constants[12])/constants[14])*constants[14] , constants[15]), constants[16] , True, 0.00000])
algebraic[69] = 0.200000/(1.00000+exp(-(states[0]-46.7000)/7.80000))
algebraic[70] = ((constants[68]*algebraic[69]*states[1])/(states[1]+constants[70]))*(states[0]-constants[69])
rates[0] = -(algebraic[46]+algebraic[48]+algebraic[50]+algebraic[54]+algebraic[57]+algebraic[58]+algebraic[68]+algebraic[60]+algebraic[61]+algebraic[62]+algebraic[63]+algebraic[64]+algebraic[65]+algebraic[66]+algebraic[70]+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[4] = 0.180640*exp(0.0357700*(states[0]+30.0000))
algebraic[14] = 0.395600*exp(-0.0623700*(states[0]+30.0000))
algebraic[5] = (0.000152000*exp(-(states[0]+13.5000)/7.00000))/(0.00670830*exp(-(states[0]+33.5000)/7.00000)+1.00000)
algebraic[15] = (0.000950000*exp((states[0]+33.5000)/7.00000))/(0.0513350*exp((states[0]+33.5000)/7.00000)+1.00000)
algebraic[6] = 1.00000/(1.00000+exp(-(states[0]+22.5000)/7.70000))
algebraic[16] = 0.493000*exp(-0.0629000*states[0])+2.05800
algebraic[7] = 1.00000/(1.00000+exp((states[0]+45.2000)/5.70000))
algebraic[17] = 270.000+1050.00/(1.00000+exp((states[0]+45.2000)/5.70000))
algebraic[8] = (4.81333e-06*(states[0]+26.5000))/(1.00000-exp(-0.128000*(states[0]+26.5000)))
algebraic[18] = 9.53333e-05*exp(-0.0380000*(states[0]+26.5000))
algebraic[19] = 0.493000*exp(-0.0629000*states[0])+2.05800
algebraic[20] = 1200.00-170.000/(1.00000+exp((states[0]+45.2000)/5.70000))
algebraic[21] = 39.3000*exp(-0.0862000*states[0])+13.1700
algebraic[10] = 0.0137330*exp(0.0381980*states[0])
algebraic[23] = 6.89000e-05*exp(-0.0417800*states[0])
algebraic[2] = 1.00000-(states[11]+states[12]+states[13]+states[14]+states[15]+states[16]+states[17])
algebraic[12] = (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[25] = 0.0500000*exp(-(states[0]+12.0000)/13.0000)
algebraic[9] = 1.00000-(states[38]+states[39]+states[37]+states[40])
algebraic[22] = 0.0223480*exp(0.0117600*states[0])
algebraic[27] = 0.0470020*exp(-0.0631000*states[0])
algebraic[28] = 0.0908210*exp(0.0233910*(states[0]+5.00000))
algebraic[32] = 0.00649700*exp(-0.0326800*(states[0]+5.00000))
algebraic[30] = (constants[45]*states[2])/(constants[46]+states[2])
algebraic[34] = 13.0000*(1.00000-exp(-(power(states[0]+14.5000, 2.00000))/100.000))
algebraic[29] = power(1.00000+(constants[18]*constants[20])/(power(constants[20]+states[3], 2.00000)), -1.00000)
algebraic[33] = constants[25]*(states[8]+states[9])*(states[3]-states[2])*states[5]
algebraic[36] = (states[4]-states[3])/constants[26]
algebraic[40] = constants[27]*(states[4]-states[1])
algebraic[42] = (constants[29]*(power(states[1], 2.00000)))/(power(constants[30], 2.00000)+power(states[1], 2.00000))
algebraic[3] = 1.00000-(states[19]+states[20]+states[21]+states[24]+states[22]+states[23]+states[25]+states[26])
algebraic[13] = 3.80200/(0.102700*exp(-(states[0]+2.50000)/17.0000)+0.200000*exp(-(states[0]+2.50000)/150.000))
algebraic[35] = 0.191700*exp(-(states[0]+2.50000)/20.3000)
algebraic[26] = 3.80200/(0.102700*exp(-(states[0]+2.50000)/15.0000)+0.230000*exp(-(states[0]+2.50000)/150.000))
algebraic[37] = 0.200000*exp(-(states[0]-2.50000)/20.3000)
algebraic[41] = 7.00000e-07*exp(-(states[0]+7.00000)/7.70000)
algebraic[43] = 0.00840000+2.00000e-05*(states[0]+7.00000)
algebraic[31] = 3.80200/(0.102700*exp(-(states[0]+2.50000)/12.0000)+0.250000*exp(-(states[0]+2.50000)/150.000))
algebraic[39] = 0.220000*exp(-(states[0]-7.50000)/20.3000)
algebraic[46] = constants[43]*states[11]*(states[0]-constants[42])
algebraic[24] = power(1.00000+(constants[17]*constants[19])/(power(constants[19]+states[2], 2.00000)), -1.00000)
algebraic[38] = (states[2]-states[1])/constants[28]
algebraic[45] = 1.00000/(0.188495*exp(-(states[0]+7.00000)/16.6000)+0.393956)
algebraic[47] = (algebraic[31]*algebraic[45]*algebraic[41])/(algebraic[39]*algebraic[43])
algebraic[49] = algebraic[45]/1000.00
algebraic[51] = algebraic[41]
algebraic[48] = (constants[47]*(power(states[1], 2.00000)))/(power(constants[48], 2.00000)+power(states[1], 2.00000))
algebraic[50] = ((((((constants[49]*1.00000)/(power(constants[50], 3.00000)+power(constants[7], 3.00000)))*1.00000)/(constants[51]+constants[8]))*1.00000)/(1.00000+constants[52]*exp(((constants[53]-1.00000)*states[0]*constants[11])/(constants[9]*constants[10]))))*(exp((constants[53]*states[0]*constants[11])/(constants[9]*constants[10]))*(power(states[18], 3.00000))*constants[8]-exp(((constants[53]-1.00000)*states[0]*constants[11])/(constants[9]*constants[10]))*(power(constants[7], 3.00000))*states[1])
algebraic[52] = ((constants[9]*constants[10])/(2.00000*constants[11]))*log(constants[8]/states[1])
algebraic[54] = constants[54]*(states[0]-algebraic[52])
algebraic[11] = power(1.00000+(constants[17]*constants[19])/(power(constants[19]+states[1], 2.00000)), -1.00000)
algebraic[44] = (constants[21]*states[1]*(constants[32]-states[7])+constants[23]*states[1]*(constants[31]-states[6]))-(constants[22]*states[7]+constants[24]*states[6])
algebraic[53] = algebraic[45]/95000.0
algebraic[55] = algebraic[41]/50.0000
algebraic[56] = ((constants[9]*constants[10])/constants[11])*log((0.900000*constants[7]+0.100000*constants[6])/(0.900000*states[18]+0.100000*states[27]))
algebraic[57] = constants[55]*states[19]*(states[0]-algebraic[56])
algebraic[58] = constants[56]*(states[0]-algebraic[56])
algebraic[67] = 1.00000/(1.00000+0.124500*exp((-0.100000*states[0]*constants[11])/(constants[9]*constants[10]))+0.0365000*constants[71]*exp((-states[0]*constants[11])/(constants[9]*constants[10])))
algebraic[68] = (((constants[65]*algebraic[67]*1.00000)/(1.00000+power(constants[66]/states[18], 1.50000)))*constants[6])/(constants[6]+constants[67])
algebraic[59] = ((constants[9]*constants[10])/constants[11])*log(constants[6]/states[27])
algebraic[60] = constants[57]*(power(states[28], 3.00000))*states[29]*(states[0]-algebraic[59])
algebraic[61] = constants[58]*states[30]*states[31]*(states[0]-algebraic[59])
algebraic[62] = (((0.293800*constants[6])/(constants[6]+210.000))*(states[0]-algebraic[59]))/(1.00000+exp(0.0896000*(states[0]-algebraic[59])))
algebraic[63] = constants[59]*(power(states[32], 2.00000))*(states[0]-algebraic[59])
algebraic[64] = constants[60]*states[33]*states[34]*(states[0]-algebraic[59])
algebraic[65] = constants[61]*states[35]*states[36]*(states[0]-algebraic[59])
algebraic[66] = constants[62]*states[37]*(states[0]-((constants[9]*constants[10])/constants[11])*log((0.980000*constants[6]+0.0200000*constants[7])/(0.980000*states[27]+0.0200000*states[18])))
algebraic[0] = custom_piecewise([greater_equal(voi , constants[12]) & less_equal(voi , constants[13]) & less_equal((voi-constants[12])-floor((voi-constants[12])/constants[14])*constants[14] , constants[15]), constants[16] , True, 0.00000])
algebraic[69] = 0.200000/(1.00000+exp(-(states[0]-46.7000)/7.80000))
algebraic[70] = ((constants[68]*algebraic[69]*states[1])/(states[1]+constants[70]))*(states[0]-constants[69])
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)