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 = 98
sizeStates = 35
sizeConstants = 103
from math import *
from numpy import *

def createLegends():
    legend_states = [""] * sizeStates
    legend_rates = [""] * sizeStates
    legend_algebraic = [""] * sizeAlgebraic
    legend_voi = ""
    legend_constants = [""] * sizeConstants
    legend_states[0] = "Ca_jsr in component Ca_dynamics (mM)"
    legend_states[1] = "Ca_sub in component Ca_dynamics (mM)"
    legend_constants[0] = "EC50_SR in component Ca_SR_release (mM)"
    legend_constants[1] = "HSR in component Ca_SR_release (dimensionless)"
    legend_states[2] = "I in component Ca_SR_release (dimensionless)"
    legend_constants[2] = "MaxSR in component Ca_SR_release (dimensionless)"
    legend_constants[3] = "MinSR in component Ca_SR_release (dimensionless)"
    legend_states[3] = "O in component Ca_SR_release (dimensionless)"
    legend_states[4] = "R in component Ca_SR_release (dimensionless)"
    legend_states[5] = "RI in component Ca_SR_release (dimensionless)"
    legend_algebraic[1] = "j_SRCarel in component Ca_SR_release (mol_per_m3_per_s)"
    legend_algebraic[0] = "kCaSR in component Ca_SR_release (dimensionless)"
    legend_constants[4] = "kiCa in component Ca_SR_release (m3_per_s_per_mol)"
    legend_algebraic[5] = "kiSRCa in component Ca_SR_release (m3_per_s_per_mol)"
    legend_constants[5] = "kim in component Ca_SR_release (hertz)"
    legend_constants[6] = "koCa in component Ca_SR_release (m6_per_s_per_mol2)"
    legend_algebraic[9] = "koSRCa in component Ca_SR_release (m6_per_s_per_mol2)"
    legend_constants[7] = "kom in component Ca_SR_release (hertz)"
    legend_constants[8] = "ks in component Ca_SR_release (hertz)"
    legend_voi = "time in component environment (second)"
    legend_constants[9] = "CM_tot in component Ca_buffering (mM)"
    legend_constants[10] = "CQ_tot in component Ca_buffering (mM)"
    legend_states[6] = "Cai in component Ca_dynamics (mM)"
    legend_constants[11] = "Mgi in component Ca_buffering (mM)"
    legend_constants[12] = "TC_tot in component Ca_buffering (mM)"
    legend_constants[13] = "TMC_tot in component Ca_buffering (mM)"
    legend_algebraic[6] = "delta_fCMi in component Ca_buffering (hertz)"
    legend_algebraic[10] = "delta_fCMs in component Ca_buffering (hertz)"
    legend_algebraic[11] = "delta_fCQ in component Ca_buffering (hertz)"
    legend_algebraic[12] = "delta_fTC in component Ca_buffering (hertz)"
    legend_algebraic[13] = "delta_fTMC in component Ca_buffering (hertz)"
    legend_algebraic[2] = "delta_fTMM in component Ca_buffering (hertz)"
    legend_states[7] = "fCMi in component Ca_buffering (dimensionless)"
    legend_states[8] = "fCMs in component Ca_buffering (dimensionless)"
    legend_states[9] = "fCQ in component Ca_buffering (dimensionless)"
    legend_states[10] = "fTC in component Ca_buffering (dimensionless)"
    legend_states[11] = "fTMC in component Ca_buffering (dimensionless)"
    legend_states[12] = "fTMM in component Ca_buffering (dimensionless)"
    legend_constants[14] = "kb_CM in component Ca_buffering (hertz)"
    legend_constants[15] = "kb_CQ in component Ca_buffering (hertz)"
    legend_constants[16] = "kb_TC in component Ca_buffering (hertz)"
    legend_constants[17] = "kb_TMC in component Ca_buffering (hertz)"
    legend_constants[18] = "kb_TMM in component Ca_buffering (hertz)"
    legend_constants[19] = "kf_CM in component Ca_buffering (m3_per_s_per_mol)"
    legend_constants[20] = "kf_CQ in component Ca_buffering (m3_per_s_per_mol)"
    legend_constants[21] = "kf_TC in component Ca_buffering (m3_per_s_per_mol)"
    legend_constants[22] = "kf_TMC in component Ca_buffering (m3_per_s_per_mol)"
    legend_constants[23] = "kf_TMM in component Ca_buffering (m3_per_s_per_mol)"
    legend_states[13] = "Ca_nsr in component Ca_dynamics (mM)"
    legend_constants[24] = "F in component Membrane (C_per_mol)"
    legend_constants[91] = "V_i in component Cell_parameters (uL)"
    legend_constants[88] = "V_jsr in component Cell_parameters (uL)"
    legend_constants[89] = "V_nsr in component Cell_parameters (uL)"
    legend_constants[90] = "V_sub in component Cell_parameters (uL)"
    legend_algebraic[64] = "i_CaT in component i_CaT (nA)"
    legend_algebraic[90] = "i_NaCa in component i_NaCa (nA)"
    legend_algebraic[20] = "i_siCa in component i_CaL (nA)"
    legend_algebraic[14] = "j_Ca_dif in component Ca_intracellular_fluxes (mol_per_m3_per_s)"
    legend_algebraic[15] = "j_tr in component Ca_intracellular_fluxes (mol_per_m3_per_s)"
    legend_algebraic[17] = "j_up in component Ca_intracellular_fluxes (mol_per_m3_per_s)"
    legend_constants[25] = "ACh in component Rate_modulation_experiments (mM)"
    legend_constants[26] = "Iso_1_uM in component Rate_modulation_experiments (dimensionless)"
    legend_constants[27] = "K_up in component Ca_intracellular_fluxes (mM)"
    legend_constants[87] = "P_up in component Ca_intracellular_fluxes (mol_per_m3_per_s)"
    legend_constants[28] = "P_up_basal in component Ca_intracellular_fluxes (mol_per_m3_per_s)"
    legend_constants[82] = "b_up in component Ca_intracellular_fluxes (dimensionless)"
    legend_constants[29] = "slope_up in component Ca_intracellular_fluxes (mM)"
    legend_constants[30] = "tau_dif_Ca in component Ca_intracellular_fluxes (second)"
    legend_constants[31] = "tau_tr in component Ca_intracellular_fluxes (second)"
    legend_constants[32] = "L_cell in component Cell_parameters (um)"
    legend_constants[33] = "L_sub in component Cell_parameters (um)"
    legend_constants[34] = "R_cell in component Cell_parameters (um)"
    legend_constants[83] = "V_cell in component Cell_parameters (uL)"
    legend_constants[35] = "V_i_part in component Cell_parameters (dimensionless)"
    legend_constants[36] = "V_jsr_part in component Cell_parameters (dimensionless)"
    legend_constants[37] = "V_nsr_part in component Cell_parameters (dimensionless)"
    legend_constants[38] = "Cao in component Ionic_values (mM)"
    legend_algebraic[16] = "E_K in component Ionic_values (mV)"
    legend_algebraic[18] = "E_Na in component Ionic_values (mV)"
    legend_states[14] = "Ki in component Ki_concentration (mM)"
    legend_constants[39] = "Ko in component Ionic_values (mM)"
    legend_states[15] = "Nai in component Nai_concentration (mM)"
    legend_constants[40] = "Nao in component Ionic_values (mM)"
    legend_constants[92] = "RTONF in component Membrane (mV)"
    legend_algebraic[69] = "i_Kr in component i_Kr (nA)"
    legend_algebraic[71] = "i_Ks in component i_Ks (nA)"
    legend_algebraic[72] = "i_Kur in component i_Kur (nA)"
    legend_algebraic[91] = "i_NaK in component i_NaK (nA)"
    legend_algebraic[92] = "i_SK in component i_SK (nA)"
    legend_algebraic[93] = "i_fK in component i_f (nA)"
    legend_algebraic[36] = "i_siK in component i_CaL (nA)"
    legend_algebraic[96] = "i_to in component i_to (nA)"
    legend_constants[41] = "C in component Membrane (uF)"
    legend_constants[42] = "R in component Membrane (mJ_per_mol_per_K)"
    legend_constants[43] = "T in component Membrane (kelvin)"
    legend_algebraic[19] = "V in component Membrane (mV)"
    legend_states[16] = "V_ode in component Membrane (mV)"
    legend_algebraic[59] = "i_CaL in component i_CaL (nA)"
    legend_algebraic[67] = "i_KACh in component i_KACh (nA)"
    legend_algebraic[76] = "i_Na in component i_Na (nA)"
    legend_algebraic[95] = "i_f in component i_f (nA)"
    legend_algebraic[97] = "i_tot in component Membrane (nA)"
    legend_algebraic[94] = "i_fNa in component i_f (nA)"
    legend_algebraic[53] = "i_siNa in component i_CaL (nA)"
    legend_constants[93] = "ACh_block in component i_CaL (dimensionless)"
    legend_constants[94] = "Iso_increase in component i_CaL (dimensionless)"
    legend_constants[44] = "P_CaL in component i_CaL (m3_A_per_mol_times_1e_minus_9)"
    legend_states[17] = "dL in component i_CaL_dL_gate (dimensionless)"
    legend_states[18] = "fCa in component i_CaL_fCa_gate (dimensionless)"
    legend_states[19] = "fL in component i_CaL_fL_gate (dimensionless)"
    legend_constants[95] = "Iso_shift_dL in component i_CaL_dL_gate (mV)"
    legend_constants[96] = "Iso_slope_dL in component i_CaL_dL_gate (dimensionless)"
    legend_constants[97] = "V_dL in component i_CaL_dL_gate (mV)"
    legend_algebraic[21] = "adVm in component i_CaL_dL_gate (mV)"
    legend_algebraic[37] = "alpha_dL in component i_CaL_dL_gate (hertz)"
    legend_algebraic[54] = "bdVm in component i_CaL_dL_gate (mV)"
    legend_algebraic[60] = "beta_dL in component i_CaL_dL_gate (hertz)"
    legend_algebraic[65] = "dL_infinity in component i_CaL_dL_gate (dimensionless)"
    legend_constants[45] = "k_dL in component i_CaL_dL_gate (mV)"
    legend_algebraic[68] = "tau_dL in component i_CaL_dL_gate (second)"
    legend_constants[46] = "Km_fCa in component i_CaL_fCa_gate (mM)"
    legend_constants[47] = "alpha_fCa in component i_CaL_fCa_gate (hertz)"
    legend_algebraic[3] = "fCa_infinity in component i_CaL_fCa_gate (dimensionless)"
    legend_algebraic[7] = "tau_fCa in component i_CaL_fCa_gate (second)"
    legend_algebraic[22] = "fL_infinity in component i_CaL_fL_gate (dimensionless)"
    legend_constants[48] = "k_fL in component i_CaL_fL_gate (mV)"
    legend_constants[49] = "shift_fL in component i_CaL_fL_gate (mV)"
    legend_algebraic[38] = "tau_fL in component i_CaL_fL_gate (second)"
    legend_constants[50] = "P_CaT in component i_CaT (m3_A_per_mol_times_1e_minus_9)"
    legend_states[20] = "dT in component i_CaT_dT_gate (dimensionless)"
    legend_states[21] = "fT in component i_CaT_fT_gate (dimensionless)"
    legend_algebraic[23] = "dT_infinity in component i_CaT_dT_gate (dimensionless)"
    legend_algebraic[39] = "tau_dT in component i_CaT_dT_gate (second)"
    legend_algebraic[24] = "fT_infinity in component i_CaT_fT_gate (dimensionless)"
    legend_constants[51] = "offset_fT in component i_CaT_fT_gate (second)"
    legend_algebraic[40] = "tau_fT in component i_CaT_fT_gate (second)"
    legend_constants[52] = "ACh_on in component i_KACh (dimensionless)"
    legend_states[22] = "a in component i_KACh_a_gate (dimensionless)"
    legend_constants[53] = "g_KACh in component i_KACh (uS)"
    legend_algebraic[41] = "a_infinity in component i_KACh_a_gate (dimensionless)"
    legend_constants[98] = "alpha_a in component i_KACh_a_gate (hertz)"
    legend_algebraic[25] = "beta_a in component i_KACh_a_gate (hertz)"
    legend_algebraic[55] = "tau_a in component i_KACh_a_gate (second)"
    legend_constants[84] = "g_Kr in component i_Kr (uS)"
    legend_states[23] = "paF in component i_Kr_pa_gate (dimensionless)"
    legend_states[24] = "paS in component i_Kr_pa_gate (dimensionless)"
    legend_states[25] = "piy in component i_Kr_pi_gate (dimensionless)"
    legend_algebraic[26] = "pa_infinity in component i_Kr_pa_gate (dimensionless)"
    legend_algebraic[42] = "tau_paF in component i_Kr_pa_gate (second)"
    legend_algebraic[43] = "tau_paS in component i_Kr_pa_gate (second)"
    legend_algebraic[27] = "pi_infinity in component i_Kr_pi_gate (dimensionless)"
    legend_algebraic[44] = "tau_pi in component i_Kr_pi_gate (second)"
    legend_algebraic[70] = "E_Ks in component i_Ks (mV)"
    legend_constants[85] = "g_Ks in component i_Ks (uS)"
    legend_constants[54] = "g_Ks_ in component i_Ks (uS)"
    legend_states[26] = "n in component i_Ks_n_gate (dimensionless)"
    legend_constants[99] = "Iso_shift in component i_Ks_n_gate (mV)"
    legend_algebraic[28] = "alpha_n in component i_Ks_n_gate (hertz)"
    legend_algebraic[45] = "beta_n in component i_Ks_n_gate (hertz)"
    legend_algebraic[56] = "n_infinity in component i_Ks_n_gate (dimensionless)"
    legend_algebraic[61] = "tau_n in component i_Ks_n_gate (second)"
    legend_constants[55] = "g_Kur in component i_Kur (uS)"
    legend_states[27] = "r_Kur in component i_Kur_rKur_gate (dimensionless)"
    legend_states[28] = "s_Kur in component i_Kur_sKur_gate (dimensionless)"
    legend_algebraic[29] = "r_Kur_infinity in component i_Kur_rKur_gate (dimensionless)"
    legend_algebraic[46] = "tau_r_Kur in component i_Kur_rKur_gate (second)"
    legend_algebraic[30] = "s_Kur_infinity in component i_Kur_sKur_gate (dimensionless)"
    legend_algebraic[47] = "tau_s_Kur in component i_Kur_sKur_gate (second)"
    legend_algebraic[73] = "E_mh in component i_Na (mV)"
    legend_constants[56] = "g_Na in component i_Na (uS)"
    legend_constants[57] = "g_Na_L in component i_Na (uS)"
    legend_states[29] = "h in component i_Na_h_gate (dimensionless)"
    legend_algebraic[74] = "i_Na_ in component i_Na (nA)"
    legend_algebraic[75] = "i_Na_L in component i_Na (nA)"
    legend_states[30] = "m in component i_Na_m_gate (dimensionless)"
    legend_constants[58] = "K1ni in component i_NaCa (mM)"
    legend_constants[59] = "K1no in component i_NaCa (mM)"
    legend_constants[60] = "K2ni in component i_NaCa (mM)"
    legend_constants[61] = "K2no in component i_NaCa (mM)"
    legend_constants[62] = "K3ni in component i_NaCa (mM)"
    legend_constants[63] = "K3no in component i_NaCa (mM)"
    legend_constants[64] = "K_NaCa in component i_NaCa (nA)"
    legend_constants[65] = "Kci in component i_NaCa (mM)"
    legend_constants[66] = "Kcni in component i_NaCa (mM)"
    legend_constants[67] = "Kco in component i_NaCa (mM)"
    legend_constants[68] = "Qci in component i_NaCa (dimensionless)"
    legend_constants[69] = "Qco in component i_NaCa (dimensionless)"
    legend_constants[70] = "Qn in component i_NaCa (dimensionless)"
    legend_algebraic[77] = "di in component i_NaCa (dimensionless)"
    legend_algebraic[78] = "do in component i_NaCa (dimensionless)"
    legend_algebraic[79] = "k12 in component i_NaCa (dimensionless)"
    legend_algebraic[80] = "k14 in component i_NaCa (dimensionless)"
    legend_algebraic[81] = "k21 in component i_NaCa (dimensionless)"
    legend_algebraic[82] = "k23 in component i_NaCa (dimensionless)"
    legend_algebraic[83] = "k32 in component i_NaCa (dimensionless)"
    legend_constants[100] = "k34 in component i_NaCa (dimensionless)"
    legend_algebraic[84] = "k41 in component i_NaCa (dimensionless)"
    legend_algebraic[85] = "k43 in component i_NaCa (dimensionless)"
    legend_algebraic[86] = "x1 in component i_NaCa (dimensionless)"
    legend_algebraic[87] = "x2 in component i_NaCa (dimensionless)"
    legend_algebraic[88] = "x3 in component i_NaCa (dimensionless)"
    legend_algebraic[89] = "x4 in component i_NaCa (dimensionless)"
    legend_constants[86] = "Iso_increase in component i_NaK (dimensionless)"
    legend_constants[71] = "Km_Kp in component i_NaK (mM)"
    legend_constants[72] = "Km_Nap in component i_NaK (mM)"
    legend_constants[73] = "i_NaK_max in component i_NaK (nA)"
    legend_algebraic[31] = "alpha_h in component i_Na_h_gate (hertz)"
    legend_algebraic[48] = "beta_h in component i_Na_h_gate (hertz)"
    legend_algebraic[57] = "h_infinity in component i_Na_h_gate (dimensionless)"
    legend_algebraic[62] = "tau_h in component i_Na_h_gate (second)"
    legend_algebraic[32] = "E0_m in component i_Na_m_gate (mV)"
    legend_algebraic[49] = "alpha_m in component i_Na_m_gate (hertz)"
    legend_algebraic[58] = "beta_m in component i_Na_m_gate (hertz)"
    legend_constants[74] = "delta_m in component i_Na_m_gate (mV)"
    legend_algebraic[63] = "m_infinity in component i_Na_m_gate (dimensionless)"
    legend_algebraic[66] = "tau_m in component i_Na_m_gate (second)"
    legend_constants[75] = "g_SK in component i_SK (uS)"
    legend_states[31] = "x in component i_SK_x_gate (dimensionless)"
    legend_constants[76] = "EC50_SK in component i_SK_x_gate (mM)"
    legend_constants[77] = "n_SK in component i_SK_x_gate (dimensionless)"
    legend_algebraic[4] = "tau_x in component i_SK_x_gate (second)"
    legend_algebraic[8] = "x_infinity in component i_SK_x_gate (dimensionless)"
    legend_constants[78] = "g_f_K in component i_f (uS)"
    legend_constants[79] = "g_f_Na in component i_f (uS)"
    legend_states[32] = "y in component i_f_y_gate (dimensionless)"
    legend_constants[101] = "ACh_shift in component i_f_y_gate (mV)"
    legend_constants[102] = "Iso_shift in component i_f_y_gate (mV)"
    legend_algebraic[33] = "tau_y in component i_f_y_gate (second)"
    legend_algebraic[50] = "y_infinity in component i_f_y_gate (dimensionless)"
    legend_constants[80] = "y_shift in component i_f_y_gate (mV)"
    legend_constants[81] = "g_to in component i_to (uS)"
    legend_states[33] = "q in component i_to_q_gate (dimensionless)"
    legend_states[34] = "r in component i_to_r_gate (dimensionless)"
    legend_algebraic[34] = "q_infinity in component i_to_q_gate (dimensionless)"
    legend_algebraic[51] = "tau_q in component i_to_q_gate (second)"
    legend_algebraic[35] = "r_infinity in component i_to_r_gate (dimensionless)"
    legend_algebraic[52] = "tau_r in component i_to_r_gate (second)"
    legend_rates[2] = "d/dt I in component Ca_SR_release (dimensionless)"
    legend_rates[3] = "d/dt O in component Ca_SR_release (dimensionless)"
    legend_rates[4] = "d/dt R in component Ca_SR_release (dimensionless)"
    legend_rates[5] = "d/dt RI in component Ca_SR_release (dimensionless)"
    legend_rates[7] = "d/dt fCMi in component Ca_buffering (dimensionless)"
    legend_rates[8] = "d/dt fCMs in component Ca_buffering (dimensionless)"
    legend_rates[9] = "d/dt fCQ in component Ca_buffering (dimensionless)"
    legend_rates[10] = "d/dt fTC in component Ca_buffering (dimensionless)"
    legend_rates[11] = "d/dt fTMC in component Ca_buffering (dimensionless)"
    legend_rates[12] = "d/dt fTMM in component Ca_buffering (dimensionless)"
    legend_rates[0] = "d/dt Ca_jsr in component Ca_dynamics (mM)"
    legend_rates[13] = "d/dt Ca_nsr in component Ca_dynamics (mM)"
    legend_rates[1] = "d/dt Ca_sub in component Ca_dynamics (mM)"
    legend_rates[6] = "d/dt Cai in component Ca_dynamics (mM)"
    legend_rates[14] = "d/dt Ki in component Ki_concentration (mM)"
    legend_rates[16] = "d/dt V_ode in component Membrane (mV)"
    legend_rates[15] = "d/dt Nai in component Nai_concentration (mM)"
    legend_rates[17] = "d/dt dL in component i_CaL_dL_gate (dimensionless)"
    legend_rates[18] = "d/dt fCa in component i_CaL_fCa_gate (dimensionless)"
    legend_rates[19] = "d/dt fL in component i_CaL_fL_gate (dimensionless)"
    legend_rates[20] = "d/dt dT in component i_CaT_dT_gate (dimensionless)"
    legend_rates[21] = "d/dt fT in component i_CaT_fT_gate (dimensionless)"
    legend_rates[22] = "d/dt a in component i_KACh_a_gate (dimensionless)"
    legend_rates[23] = "d/dt paF in component i_Kr_pa_gate (dimensionless)"
    legend_rates[24] = "d/dt paS in component i_Kr_pa_gate (dimensionless)"
    legend_rates[25] = "d/dt piy in component i_Kr_pi_gate (dimensionless)"
    legend_rates[26] = "d/dt n in component i_Ks_n_gate (dimensionless)"
    legend_rates[27] = "d/dt r_Kur in component i_Kur_rKur_gate (dimensionless)"
    legend_rates[28] = "d/dt s_Kur in component i_Kur_sKur_gate (dimensionless)"
    legend_rates[29] = "d/dt h in component i_Na_h_gate (dimensionless)"
    legend_rates[30] = "d/dt m in component i_Na_m_gate (dimensionless)"
    legend_rates[31] = "d/dt x in component i_SK_x_gate (dimensionless)"
    legend_rates[32] = "d/dt y in component i_f_y_gate (dimensionless)"
    legend_rates[33] = "d/dt q in component i_to_q_gate (dimensionless)"
    legend_rates[34] = "d/dt r in component i_to_r_gate (dimensionless)"
    return (legend_states, legend_algebraic, legend_voi, legend_constants)

def initConsts():
    constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
    states[0] = 0.5869378
    states[1] = 8.59587700000000057e-05
    constants[0] = 0.45
    constants[1] = 2.5
    states[2] = 1.88065699999999987e-09
    constants[2] = 15.0
    constants[3] = 1.0
    states[3] = 1.75000900000000001e-08
    states[4] = 0.9004965
    states[5] = 9.67726599999999965e-02
    constants[4] = 500.0
    constants[5] = 5.0
    constants[6] = 10000.0
    constants[7] = 660.0
    constants[8] = 1.48041085099999994e+08
    constants[9] = 0.045
    constants[10] = 10.0
    states[6] = 1.26704600000000013e-04
    constants[11] = 2.5
    constants[12] = 0.031
    constants[13] = 0.062
    states[7] = 0.2776014
    states[8] = 0.206602
    states[9] = 0.1878388
    states[10] = 0.0246504
    states[11] = 0.329845
    states[12] = 0.5920168
    constants[14] = 542.0
    constants[15] = 445.0
    constants[16] = 446.0
    constants[17] = 7.51
    constants[18] = 751.0
    constants[19] = 1641986.0
    constants[20] = 175.4
    constants[21] = 88800.0
    constants[22] = 227700.0
    constants[23] = 2277.0
    states[13] = 0.669119
    constants[24] = 9.64853414999999950e+04
    constants[25] = 0.0
    constants[26] = 0.0
    constants[27] = 2.86113000000000003e-04
    constants[28] = 5.0
    constants[29] = 5e-05
    constants[30] = 5.469e-05
    constants[31] = 0.04
    constants[32] = 67.0
    constants[33] = 0.02
    constants[34] = 3.9
    constants[35] = 0.46
    constants[36] = 0.0012
    constants[37] = 0.0116
    constants[38] = 1.8
    states[14] = 139.1382
    constants[39] = 5.4
    states[15] = 6.084085
    constants[40] = 140.0
    constants[41] = 5.7e-05
    constants[42] = 8314.472
    constants[43] = 310.0
    states[16] = -4.61898999999999990e-02
    constants[44] = 0.5046812
    states[17] = 1.40229900000000008e-03
    states[18] = 0.7820265
    states[19] = 0.9951828
    constants[45] = 5.854
    constants[46] = 0.000338
    constants[47] = 0.0075
    constants[48] = 5.1737
    constants[49] = 0.0
    constants[50] = 0.04132
    states[20] = 0.192108
    states[21] = 3.74674099999999996e-02
    constants[51] = 0.0
    constants[52] = 0.0
    states[22] = 2.81845900000000016e-03
    constants[53] = 0.00345
    states[23] = 8.34449800000000054e-03
    states[24] = 0.3450797
    states[25] = 0.7369347
    constants[54] = 8.63714000000000018e-04
    states[26] = 0.1057902
    constants[55] = 7.062e-05
    states[27] = 9.05902400000000059e-03
    states[28] = 0.8655135
    constants[56] = 0.0223
    constants[57] = 0.0
    states[29] = 5.20200199999999984e-03
    states[30] = 0.3777047
    constants[58] = 395.3
    constants[59] = 1628.0
    constants[60] = 2.289
    constants[61] = 561.4
    constants[62] = 26.44
    constants[63] = 4.663
    constants[64] = 3.343
    constants[65] = 0.0207
    constants[66] = 26.44
    constants[67] = 3.663
    constants[68] = 0.1369
    constants[69] = 0.0
    constants[70] = 0.4315
    constants[71] = 1.4
    constants[72] = 14.0
    constants[73] = 0.137171
    constants[74] = 1e-05
    constants[75] = 0.000165
    states[31] = 6.77693899999999988e-02
    constants[76] = 0.7
    constants[77] = 2.2
    constants[78] = 0.0117
    constants[79] = 0.00696
    states[32] = 0.0103301
    constants[80] = 0.0
    constants[81] = 1.67347999999999998e-03
    states[33] = 0.4735559
    states[34] = 1.24570600000000007e-02
    constants[82] = custom_piecewise([greater(constants[26] , 0.00000), -0.250000 , greater(constants[25] , 0.00000), (0.700000*constants[25])/(9.00000e-05+constants[25]) , True, 0.00000])
    constants[83] = ((1.00000e-09*3.14159)*(power(constants[34], 2.00000)))*constants[32]
    constants[84] = 0.00498910*(power(constants[39]/5.40000, 1.0/2))
    constants[85] = custom_piecewise([greater(constants[26] , 0.00000), 1.20000*constants[54] , True, constants[54]])
    constants[86] = custom_piecewise([greater(constants[26] , 0.00000), 1.20000 , True, 1.00000])
    constants[87] = constants[28]*(1.00000-constants[82])
    constants[88] = constants[36]*constants[83]
    constants[89] = constants[37]*constants[83]
    constants[90] = ((((1.00000e-09*2.00000)*3.14159)*constants[33])*(constants[34]-constants[33]/2.00000))*constants[32]
    constants[91] = constants[35]*constants[83]-constants[90]
    constants[92] = (constants[42]*constants[43])/constants[24]
    constants[93] = (0.310000*constants[25])/(constants[25]+9.00000e-05)
    constants[94] = custom_piecewise([greater(constants[26] , 0.00000), 1.23000 , True, 1.00000])
    constants[95] = custom_piecewise([greater(constants[26] , 0.00000), -8.00000 , True, 0.00000])
    constants[96] = custom_piecewise([greater(constants[26] , 0.00000), -27.0000 , True, 0.00000])
    constants[97] = -7.77130
    constants[98] = (3.59880-0.0256410)/(1.00000+1.21550e-06/(power(1.00000*constants[25], 1.69510)))+0.0256410
    constants[99] = custom_piecewise([greater(constants[26] , 0.00000), -14.0000 , True, 0.00000])
    constants[100] = constants[40]/(constants[63]+constants[40])
    constants[101] = custom_piecewise([greater(constants[25] , 0.00000), -1.00000-(9.89800*(power(1.00000*constants[25], 0.618000)))/(power(1.00000*constants[25], 0.618000)+0.00122423) , True, 0.00000])
    constants[102] = custom_piecewise([greater(constants[26] , 0.00000), 7.50000 , True, 0.00000])
    return (states, constants)

def computeRates(voi, states, constants):
    rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
    algebraic[2] = (constants[23]*constants[11])*(1.00000-(states[11]+states[12]))-constants[18]*states[12]
    rates[12] = algebraic[2]
    algebraic[6] = (constants[19]*states[6])*(1.00000-states[7])-constants[14]*states[7]
    rates[7] = algebraic[6]
    algebraic[3] = constants[46]/(constants[46]+states[1])
    algebraic[7] = (0.00100000*algebraic[3])/constants[47]
    rates[18] = (algebraic[3]-states[18])/algebraic[7]
    algebraic[4] = 0.00100000/(0.0470000*((1000.00*states[1])/1.00000)+1.00000/76.0000)
    algebraic[8] = (0.810000*(power(states[1]/1.00000, constants[77])))/(power(states[1]/1.00000, constants[77])+power(constants[76]/1.00000, constants[77]))
    rates[31] = (algebraic[8]-states[31])/algebraic[4]
    algebraic[0] = constants[2]-(constants[2]-constants[3])/(1.00000+power(constants[0]/states[0], constants[1]))
    algebraic[5] = constants[4]*algebraic[0]
    algebraic[9] = constants[6]/algebraic[0]
    rates[2] = ((algebraic[5]*states[1])*states[3]-constants[5]*states[2])-(constants[7]*states[2]-(algebraic[9]*(power(states[1], 2.00000)))*states[5])
    rates[3] = ((algebraic[9]*(power(states[1], 2.00000)))*states[4]-constants[7]*states[3])-((algebraic[5]*states[1])*states[3]-constants[5]*states[2])
    rates[4] = (constants[5]*states[5]-(algebraic[5]*states[1])*states[4])-((algebraic[9]*(power(states[1], 2.00000)))*states[4]-constants[7]*states[3])
    rates[5] = (constants[7]*states[2]-(algebraic[9]*(power(states[1], 2.00000)))*states[5])-(constants[5]*states[5]-(algebraic[5]*states[1])*states[4])
    algebraic[10] = (constants[19]*states[1])*(1.00000-states[8])-constants[14]*states[8]
    rates[8] = algebraic[10]
    algebraic[11] = (constants[20]*states[0])*(1.00000-states[9])-constants[15]*states[9]
    rates[9] = algebraic[11]
    algebraic[12] = (constants[21]*states[6])*(1.00000-states[10])-constants[16]*states[10]
    rates[10] = algebraic[12]
    algebraic[13] = (constants[22]*states[6])*(1.00000-(states[11]+states[12]))-constants[17]*states[11]
    rates[11] = algebraic[13]
    algebraic[1] = (constants[8]*states[3])*(states[0]-states[1])
    algebraic[15] = (states[13]-states[0])/constants[31]
    rates[0] = algebraic[15]-(algebraic[1]+constants[10]*algebraic[11])
    algebraic[17] = constants[87]/(1.00000+exp((-states[6]+constants[27])/constants[29]))
    rates[13] = algebraic[17]-(algebraic[15]*constants[88])/constants[89]
    algebraic[14] = (states[1]-states[6])/constants[30]
    rates[6] = (1.00000*(algebraic[14]*constants[90]-algebraic[17]*constants[89]))/constants[91]-((constants[9]*algebraic[6]+constants[12]*algebraic[12])+constants[13]*algebraic[13])
    algebraic[19] = states[16]
    algebraic[22] = 1.00000/(1.00000+exp(((algebraic[19]+12.1847)+constants[49])/(5.30000+constants[48])))
    algebraic[38] = 0.00100000*(44.3000+230.000*exp(-(power((algebraic[19]+36.0000)/10.0000, 2.00000))))
    rates[19] = (algebraic[22]-states[19])/algebraic[38]
    algebraic[23] = 1.00000/(1.00000+exp(-(algebraic[19]+38.3000)/5.50000))
    algebraic[39] = 0.00100000/(1.06800*exp((algebraic[19]+38.3000)/30.0000)+1.06800*exp(-(algebraic[19]+38.3000)/30.0000))
    rates[20] = (algebraic[23]-states[20])/algebraic[39]
    algebraic[24] = 1.00000/(1.00000+exp((algebraic[19]+58.7000)/3.80000))
    algebraic[40] = 1.00000/(16.6700*exp(-(algebraic[19]+75.0000)/83.3000)+16.6700*exp((algebraic[19]+75.0000)/15.3800))+constants[51]
    rates[21] = (algebraic[24]-states[21])/algebraic[40]
    algebraic[26] = 1.00000/(1.00000+exp(-(algebraic[19]+10.0144)/7.66070))
    algebraic[42] = 1.00000/(30.0000*exp(algebraic[19]/10.0000)+exp(-algebraic[19]/12.0000))
    rates[23] = (algebraic[26]-states[23])/algebraic[42]
    algebraic[43] = 0.846554/(4.20000*exp(algebraic[19]/17.0000)+0.150000*exp(-algebraic[19]/21.6000))
    rates[24] = (algebraic[26]-states[24])/algebraic[43]
    algebraic[27] = 1.00000/(1.00000+exp((algebraic[19]+28.6000)/17.1000))
    algebraic[44] = 1.00000/(100.000*exp(-algebraic[19]/54.6450)+656.000*exp(algebraic[19]/106.157))
    rates[25] = (algebraic[27]-states[25])/algebraic[44]
    algebraic[29] = 1.00000/(1.00000+exp((algebraic[19]+6.00000)/-8.60000))
    algebraic[46] = 0.00900000/(1.00000+exp((algebraic[19]+5.00000)/12.0000))+0.000500000
    rates[27] = (algebraic[29]-states[27])/algebraic[46]
    algebraic[30] = 1.00000/(1.00000+exp((algebraic[19]+7.50000)/10.0000))
    algebraic[47] = 0.590000/(1.00000+exp((algebraic[19]+60.0000)/10.0000))+3.05000
    rates[28] = (algebraic[30]-states[28])/algebraic[47]
    algebraic[33] = 1.00000/((0.360000*(((algebraic[19]+148.800)-constants[101])-constants[102]))/(exp(0.0660000*(((algebraic[19]+148.800)-constants[101])-constants[102]))-1.00000)+(0.100000*(((algebraic[19]+87.3000)-constants[101])-constants[102]))/(1.00000-exp(-0.200000*(((algebraic[19]+87.3000)-constants[101])-constants[102]))))-0.0540000
    algebraic[50] = custom_piecewise([less(algebraic[19] , -(((80.0000-constants[101])-constants[102])-constants[80])), 0.0132900+0.999210/(1.00000+exp(((((algebraic[19]+97.1340)-constants[101])-constants[102])-constants[80])/8.17520)) , True, 0.000250100*exp(-(((algebraic[19]-constants[101])-constants[102])-constants[80])/12.8610)])
    rates[32] = (algebraic[50]-states[32])/algebraic[33]
    algebraic[34] = 1.00000/(1.00000+exp((algebraic[19]+49.0000)/13.0000))
    algebraic[51] = (0.00100000*0.600000)*(65.1700/(0.570000*exp(-0.0800000*(algebraic[19]+44.0000))+0.0650000*exp(0.100000*(algebraic[19]+45.9300)))+10.1000)
    rates[33] = (algebraic[34]-states[33])/algebraic[51]
    algebraic[35] = 1.00000/(1.00000+exp(-(algebraic[19]-19.3000)/15.0000))
    algebraic[52] = ((0.00100000*0.660000)*1.40000)*(15.5900/(1.03700*exp(0.0900000*(algebraic[19]+30.6100))+0.369000*exp(-0.120000*(algebraic[19]+23.8400)))+2.98000)
    rates[34] = (algebraic[35]-states[34])/algebraic[52]
    algebraic[25] = 10.0000*exp(0.0133000*(algebraic[19]+40.0000))
    algebraic[41] = constants[98]/(constants[98]+algebraic[25])
    algebraic[55] = 1.00000/(constants[98]+algebraic[25])
    rates[22] = (algebraic[41]-states[22])/algebraic[55]
    algebraic[56] = power(1.00000/(1.00000+exp(-((algebraic[19]+0.638300)-constants[99])/10.7071)), 1.0/2)
    algebraic[28] = 28.0000/(1.00000+exp(-((algebraic[19]-40.0000)-constants[99])/3.00000))
    algebraic[45] = 1.00000*exp(-((algebraic[19]-constants[99])-5.00000)/25.0000)
    algebraic[61] = 1.00000/(algebraic[28]+algebraic[45])
    rates[26] = (algebraic[56]-states[26])/algebraic[61]
    algebraic[57] = 1.00000/(1.00000+exp((algebraic[19]+69.8040)/4.45650))
    algebraic[31] = 20.0000*exp(-0.125000*(algebraic[19]+75.0000))
    algebraic[48] = 2000.00/(320.000*exp(-0.100000*(algebraic[19]+75.0000))+1.00000)
    algebraic[62] = 1.00000/(algebraic[31]+algebraic[48])
    rates[29] = (algebraic[57]-states[29])/algebraic[62]
    algebraic[63] = 1.00000/(1.00000+exp(-(algebraic[19]+42.0504)/8.31060))
    algebraic[32] = algebraic[19]+41.0000
    algebraic[49] = custom_piecewise([less(fabs(algebraic[32]) , constants[74]), 2000.00 , True, (200.000*algebraic[32])/(1.00000-exp(-0.100000*algebraic[32]))])
    algebraic[58] = 8000.00*exp(-0.0560000*(algebraic[19]+66.0000))
    algebraic[66] = 1.00000/(algebraic[49]+algebraic[58])
    rates[30] = (algebraic[63]-states[30])/algebraic[66]
    algebraic[65] = 1.00000/(1.00000+exp(-((algebraic[19]-constants[97])-constants[95])/(constants[45]*(1.00000+constants[96]/100.000))))
    algebraic[21] = custom_piecewise([equal(algebraic[19] , -41.8000), -41.8000 , equal(algebraic[19] , 0.00000), 0.00000 , equal(algebraic[19] , -6.80000), -6.80001 , True, algebraic[19]])
    algebraic[37] = (-0.0283900*(algebraic[21]+41.8000))/(exp(-(algebraic[21]+41.8000)/2.50000)-1.00000)-(0.0849000*(algebraic[21]+6.80000))/(exp(-(algebraic[21]+6.80000)/4.80000)-1.00000)
    algebraic[54] = custom_piecewise([equal(algebraic[19] , -1.80000), -1.80001 , True, algebraic[19]])
    algebraic[60] = (0.0114300*(algebraic[54]+1.80000))/(exp((algebraic[54]+1.80000)/2.50000)-1.00000)
    algebraic[68] = 0.00100000/(algebraic[37]+algebraic[60])
    rates[17] = (algebraic[65]-states[17])/algebraic[68]
    algebraic[64] = (((((2.00000*constants[50])*algebraic[19])/(constants[92]*(1.00000-exp(((-1.00000*algebraic[19])*2.00000)/constants[92]))))*(states[1]-constants[38]*exp((-2.00000*algebraic[19])/constants[92])))*states[20])*states[21]
    algebraic[77] = (1.00000+(states[1]/constants[65])*((1.00000+exp((-constants[68]*algebraic[19])/constants[92]))+states[15]/constants[66]))+(states[15]/constants[58])*(1.00000+(states[15]/constants[60])*(1.00000+states[15]/constants[62]))
    algebraic[79] = ((states[1]/constants[65])*exp((-constants[68]*algebraic[19])/constants[92]))/algebraic[77]
    algebraic[78] = (1.00000+(constants[38]/constants[67])*(1.00000+exp((constants[69]*algebraic[19])/constants[92])))+(constants[40]/constants[59])*(1.00000+(constants[40]/constants[61])*(1.00000+constants[40]/constants[63]))
    algebraic[81] = ((constants[38]/constants[67])*exp((constants[69]*algebraic[19])/constants[92]))/algebraic[78]
    algebraic[82] = (((((constants[40]/constants[59])*constants[40])/constants[61])*(1.00000+constants[40]/constants[63]))*exp((-constants[70]*algebraic[19])/(2.00000*constants[92])))/algebraic[78]
    algebraic[83] = exp((constants[70]*algebraic[19])/(2.00000*constants[92]))
    algebraic[84] = exp((-constants[70]*algebraic[19])/(2.00000*constants[92]))
    algebraic[85] = states[15]/(constants[62]+states[15])
    algebraic[86] = (algebraic[84]*constants[100])*(algebraic[82]+algebraic[81])+(algebraic[81]*algebraic[83])*(algebraic[85]+algebraic[84])
    algebraic[80] = (((((states[15]/constants[58])*states[15])/constants[60])*(1.00000+states[15]/constants[62]))*exp((constants[70]*algebraic[19])/(2.00000*constants[92])))/algebraic[77]
    algebraic[87] = (algebraic[83]*algebraic[85])*(algebraic[80]+algebraic[79])+(algebraic[84]*algebraic[79])*(constants[100]+algebraic[83])
    algebraic[88] = (algebraic[80]*algebraic[85])*(algebraic[82]+algebraic[81])+(algebraic[79]*algebraic[82])*(algebraic[85]+algebraic[84])
    algebraic[89] = (algebraic[82]*constants[100])*(algebraic[80]+algebraic[79])+(algebraic[80]*algebraic[81])*(constants[100]+algebraic[83])
    algebraic[90] = (constants[64]*(algebraic[87]*algebraic[81]-algebraic[86]*algebraic[79]))/(((algebraic[86]+algebraic[87])+algebraic[88])+algebraic[89])
    algebraic[20] = ((((((2.00000*constants[44])*(algebraic[19]-0.00000))/(constants[92]*(1.00000-exp(((-1.00000*(algebraic[19]-0.00000))*2.00000)/constants[92]))))*(states[1]-constants[38]*exp((-2.00000*(algebraic[19]-0.00000))/constants[92])))*states[17])*states[19])*states[18]
    rates[1] = (algebraic[1]*constants[88])/constants[90]-((((algebraic[20]+algebraic[64])-2.00000*algebraic[90])/((2.00000*constants[24])*constants[90])+algebraic[14])+constants[9]*algebraic[10])
    algebraic[18] = constants[92]*log(constants[40]/states[15])
    algebraic[91] = (((constants[86]*constants[73])*(power(1.00000+power(constants[71]/constants[39], 1.20000), -1.00000)))*(power(1.00000+power(constants[72]/states[15], 1.30000), -1.00000)))*(power(1.00000+exp(-((algebraic[19]-algebraic[18])+110.000)/20.0000), -1.00000))
    algebraic[73] = constants[92]*log((constants[40]+0.120000*constants[39])/(states[15]+0.120000*states[14]))
    algebraic[74] = ((constants[56]*(power(states[30], 3.00000)))*states[29])*(algebraic[19]-algebraic[73])
    algebraic[75] = (constants[57]*(power(states[30], 3.00000)))*(algebraic[19]-algebraic[73])
    algebraic[76] = algebraic[74]+algebraic[75]
    algebraic[94] = (states[32]*constants[79])*(algebraic[19]-algebraic[18])
    algebraic[53] = ((((((1.85000e-05*constants[44])*(algebraic[19]-0.00000))/(constants[92]*(1.00000-exp((-1.00000*(algebraic[19]-0.00000))/constants[92]))))*(states[15]-constants[40]*exp((-1.00000*(algebraic[19]-0.00000))/constants[92])))*states[17])*states[19])*states[18]
    rates[15] = (-1.00000*((((algebraic[76]+algebraic[94])+algebraic[53])+3.00000*algebraic[91])+3.00000*algebraic[90]))/((1.00000*(constants[91]+constants[90]))*constants[24])
    algebraic[16] = constants[92]*log(constants[39]/states[14])
    algebraic[69] = ((constants[84]*(algebraic[19]-algebraic[16]))*(0.900000*states[23]+0.100000*states[24]))*states[25]
    algebraic[70] = constants[92]*log((constants[39]+0.120000*constants[40])/(states[14]+0.120000*states[15]))
    algebraic[71] = (constants[85]*(algebraic[19]-algebraic[70]))*(power(states[26], 2.00000))
    algebraic[72] = ((constants[55]*states[27])*states[28])*(algebraic[19]-algebraic[16])
    algebraic[92] = (constants[75]*(algebraic[19]-algebraic[16]))*states[31]
    algebraic[93] = (states[32]*constants[78])*(algebraic[19]-algebraic[16])
    algebraic[36] = ((((((0.000365000*constants[44])*(algebraic[19]-0.00000))/(constants[92]*(1.00000-exp((-1.00000*(algebraic[19]-0.00000))/constants[92]))))*(states[14]-constants[39]*exp((-1.00000*(algebraic[19]-0.00000))/constants[92])))*states[17])*states[19])*states[18]
    algebraic[96] = ((constants[81]*(algebraic[19]-algebraic[16]))*states[33])*states[34]
    rates[14] = (-1.00000*(((((((algebraic[72]+algebraic[96])+algebraic[69])+algebraic[71])+algebraic[93])+algebraic[36])+algebraic[92])-2.00000*algebraic[91]))/((1.00000*(constants[91]+constants[90]))*constants[24])
    algebraic[59] = ((((algebraic[20]+algebraic[36])+algebraic[53])*(1.00000-constants[93]))*1.00000)*constants[94]
    algebraic[67] = custom_piecewise([greater(constants[25] , 0.00000), (((constants[52]*constants[53])*(algebraic[19]-algebraic[16]))*(1.00000+exp((algebraic[19]+20.0000)/20.0000)))*states[22] , True, 0.00000])
    algebraic[95] = algebraic[94]+algebraic[93]
    algebraic[97] = ((((((((((algebraic[95]+algebraic[69])+algebraic[71])+algebraic[96])+algebraic[91])+algebraic[90])+algebraic[76])+algebraic[59])+algebraic[64])+algebraic[67])+algebraic[72])+algebraic[92]
    rates[16] = -algebraic[97]/constants[41]
    return(rates)

def computeAlgebraic(constants, states, voi):
    algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
    states = array(states)
    voi = array(voi)
    algebraic[2] = (constants[23]*constants[11])*(1.00000-(states[11]+states[12]))-constants[18]*states[12]
    algebraic[6] = (constants[19]*states[6])*(1.00000-states[7])-constants[14]*states[7]
    algebraic[3] = constants[46]/(constants[46]+states[1])
    algebraic[7] = (0.00100000*algebraic[3])/constants[47]
    algebraic[4] = 0.00100000/(0.0470000*((1000.00*states[1])/1.00000)+1.00000/76.0000)
    algebraic[8] = (0.810000*(power(states[1]/1.00000, constants[77])))/(power(states[1]/1.00000, constants[77])+power(constants[76]/1.00000, constants[77]))
    algebraic[0] = constants[2]-(constants[2]-constants[3])/(1.00000+power(constants[0]/states[0], constants[1]))
    algebraic[5] = constants[4]*algebraic[0]
    algebraic[9] = constants[6]/algebraic[0]
    algebraic[10] = (constants[19]*states[1])*(1.00000-states[8])-constants[14]*states[8]
    algebraic[11] = (constants[20]*states[0])*(1.00000-states[9])-constants[15]*states[9]
    algebraic[12] = (constants[21]*states[6])*(1.00000-states[10])-constants[16]*states[10]
    algebraic[13] = (constants[22]*states[6])*(1.00000-(states[11]+states[12]))-constants[17]*states[11]
    algebraic[1] = (constants[8]*states[3])*(states[0]-states[1])
    algebraic[15] = (states[13]-states[0])/constants[31]
    algebraic[17] = constants[87]/(1.00000+exp((-states[6]+constants[27])/constants[29]))
    algebraic[14] = (states[1]-states[6])/constants[30]
    algebraic[19] = states[16]
    algebraic[22] = 1.00000/(1.00000+exp(((algebraic[19]+12.1847)+constants[49])/(5.30000+constants[48])))
    algebraic[38] = 0.00100000*(44.3000+230.000*exp(-(power((algebraic[19]+36.0000)/10.0000, 2.00000))))
    algebraic[23] = 1.00000/(1.00000+exp(-(algebraic[19]+38.3000)/5.50000))
    algebraic[39] = 0.00100000/(1.06800*exp((algebraic[19]+38.3000)/30.0000)+1.06800*exp(-(algebraic[19]+38.3000)/30.0000))
    algebraic[24] = 1.00000/(1.00000+exp((algebraic[19]+58.7000)/3.80000))
    algebraic[40] = 1.00000/(16.6700*exp(-(algebraic[19]+75.0000)/83.3000)+16.6700*exp((algebraic[19]+75.0000)/15.3800))+constants[51]
    algebraic[26] = 1.00000/(1.00000+exp(-(algebraic[19]+10.0144)/7.66070))
    algebraic[42] = 1.00000/(30.0000*exp(algebraic[19]/10.0000)+exp(-algebraic[19]/12.0000))
    algebraic[43] = 0.846554/(4.20000*exp(algebraic[19]/17.0000)+0.150000*exp(-algebraic[19]/21.6000))
    algebraic[27] = 1.00000/(1.00000+exp((algebraic[19]+28.6000)/17.1000))
    algebraic[44] = 1.00000/(100.000*exp(-algebraic[19]/54.6450)+656.000*exp(algebraic[19]/106.157))
    algebraic[29] = 1.00000/(1.00000+exp((algebraic[19]+6.00000)/-8.60000))
    algebraic[46] = 0.00900000/(1.00000+exp((algebraic[19]+5.00000)/12.0000))+0.000500000
    algebraic[30] = 1.00000/(1.00000+exp((algebraic[19]+7.50000)/10.0000))
    algebraic[47] = 0.590000/(1.00000+exp((algebraic[19]+60.0000)/10.0000))+3.05000
    algebraic[33] = 1.00000/((0.360000*(((algebraic[19]+148.800)-constants[101])-constants[102]))/(exp(0.0660000*(((algebraic[19]+148.800)-constants[101])-constants[102]))-1.00000)+(0.100000*(((algebraic[19]+87.3000)-constants[101])-constants[102]))/(1.00000-exp(-0.200000*(((algebraic[19]+87.3000)-constants[101])-constants[102]))))-0.0540000
    algebraic[50] = custom_piecewise([less(algebraic[19] , -(((80.0000-constants[101])-constants[102])-constants[80])), 0.0132900+0.999210/(1.00000+exp(((((algebraic[19]+97.1340)-constants[101])-constants[102])-constants[80])/8.17520)) , True, 0.000250100*exp(-(((algebraic[19]-constants[101])-constants[102])-constants[80])/12.8610)])
    algebraic[34] = 1.00000/(1.00000+exp((algebraic[19]+49.0000)/13.0000))
    algebraic[51] = (0.00100000*0.600000)*(65.1700/(0.570000*exp(-0.0800000*(algebraic[19]+44.0000))+0.0650000*exp(0.100000*(algebraic[19]+45.9300)))+10.1000)
    algebraic[35] = 1.00000/(1.00000+exp(-(algebraic[19]-19.3000)/15.0000))
    algebraic[52] = ((0.00100000*0.660000)*1.40000)*(15.5900/(1.03700*exp(0.0900000*(algebraic[19]+30.6100))+0.369000*exp(-0.120000*(algebraic[19]+23.8400)))+2.98000)
    algebraic[25] = 10.0000*exp(0.0133000*(algebraic[19]+40.0000))
    algebraic[41] = constants[98]/(constants[98]+algebraic[25])
    algebraic[55] = 1.00000/(constants[98]+algebraic[25])
    algebraic[56] = power(1.00000/(1.00000+exp(-((algebraic[19]+0.638300)-constants[99])/10.7071)), 1.0/2)
    algebraic[28] = 28.0000/(1.00000+exp(-((algebraic[19]-40.0000)-constants[99])/3.00000))
    algebraic[45] = 1.00000*exp(-((algebraic[19]-constants[99])-5.00000)/25.0000)
    algebraic[61] = 1.00000/(algebraic[28]+algebraic[45])
    algebraic[57] = 1.00000/(1.00000+exp((algebraic[19]+69.8040)/4.45650))
    algebraic[31] = 20.0000*exp(-0.125000*(algebraic[19]+75.0000))
    algebraic[48] = 2000.00/(320.000*exp(-0.100000*(algebraic[19]+75.0000))+1.00000)
    algebraic[62] = 1.00000/(algebraic[31]+algebraic[48])
    algebraic[63] = 1.00000/(1.00000+exp(-(algebraic[19]+42.0504)/8.31060))
    algebraic[32] = algebraic[19]+41.0000
    algebraic[49] = custom_piecewise([less(fabs(algebraic[32]) , constants[74]), 2000.00 , True, (200.000*algebraic[32])/(1.00000-exp(-0.100000*algebraic[32]))])
    algebraic[58] = 8000.00*exp(-0.0560000*(algebraic[19]+66.0000))
    algebraic[66] = 1.00000/(algebraic[49]+algebraic[58])
    algebraic[65] = 1.00000/(1.00000+exp(-((algebraic[19]-constants[97])-constants[95])/(constants[45]*(1.00000+constants[96]/100.000))))
    algebraic[21] = custom_piecewise([equal(algebraic[19] , -41.8000), -41.8000 , equal(algebraic[19] , 0.00000), 0.00000 , equal(algebraic[19] , -6.80000), -6.80001 , True, algebraic[19]])
    algebraic[37] = (-0.0283900*(algebraic[21]+41.8000))/(exp(-(algebraic[21]+41.8000)/2.50000)-1.00000)-(0.0849000*(algebraic[21]+6.80000))/(exp(-(algebraic[21]+6.80000)/4.80000)-1.00000)
    algebraic[54] = custom_piecewise([equal(algebraic[19] , -1.80000), -1.80001 , True, algebraic[19]])
    algebraic[60] = (0.0114300*(algebraic[54]+1.80000))/(exp((algebraic[54]+1.80000)/2.50000)-1.00000)
    algebraic[68] = 0.00100000/(algebraic[37]+algebraic[60])
    algebraic[64] = (((((2.00000*constants[50])*algebraic[19])/(constants[92]*(1.00000-exp(((-1.00000*algebraic[19])*2.00000)/constants[92]))))*(states[1]-constants[38]*exp((-2.00000*algebraic[19])/constants[92])))*states[20])*states[21]
    algebraic[77] = (1.00000+(states[1]/constants[65])*((1.00000+exp((-constants[68]*algebraic[19])/constants[92]))+states[15]/constants[66]))+(states[15]/constants[58])*(1.00000+(states[15]/constants[60])*(1.00000+states[15]/constants[62]))
    algebraic[79] = ((states[1]/constants[65])*exp((-constants[68]*algebraic[19])/constants[92]))/algebraic[77]
    algebraic[78] = (1.00000+(constants[38]/constants[67])*(1.00000+exp((constants[69]*algebraic[19])/constants[92])))+(constants[40]/constants[59])*(1.00000+(constants[40]/constants[61])*(1.00000+constants[40]/constants[63]))
    algebraic[81] = ((constants[38]/constants[67])*exp((constants[69]*algebraic[19])/constants[92]))/algebraic[78]
    algebraic[82] = (((((constants[40]/constants[59])*constants[40])/constants[61])*(1.00000+constants[40]/constants[63]))*exp((-constants[70]*algebraic[19])/(2.00000*constants[92])))/algebraic[78]
    algebraic[83] = exp((constants[70]*algebraic[19])/(2.00000*constants[92]))
    algebraic[84] = exp((-constants[70]*algebraic[19])/(2.00000*constants[92]))
    algebraic[85] = states[15]/(constants[62]+states[15])
    algebraic[86] = (algebraic[84]*constants[100])*(algebraic[82]+algebraic[81])+(algebraic[81]*algebraic[83])*(algebraic[85]+algebraic[84])
    algebraic[80] = (((((states[15]/constants[58])*states[15])/constants[60])*(1.00000+states[15]/constants[62]))*exp((constants[70]*algebraic[19])/(2.00000*constants[92])))/algebraic[77]
    algebraic[87] = (algebraic[83]*algebraic[85])*(algebraic[80]+algebraic[79])+(algebraic[84]*algebraic[79])*(constants[100]+algebraic[83])
    algebraic[88] = (algebraic[80]*algebraic[85])*(algebraic[82]+algebraic[81])+(algebraic[79]*algebraic[82])*(algebraic[85]+algebraic[84])
    algebraic[89] = (algebraic[82]*constants[100])*(algebraic[80]+algebraic[79])+(algebraic[80]*algebraic[81])*(constants[100]+algebraic[83])
    algebraic[90] = (constants[64]*(algebraic[87]*algebraic[81]-algebraic[86]*algebraic[79]))/(((algebraic[86]+algebraic[87])+algebraic[88])+algebraic[89])
    algebraic[20] = ((((((2.00000*constants[44])*(algebraic[19]-0.00000))/(constants[92]*(1.00000-exp(((-1.00000*(algebraic[19]-0.00000))*2.00000)/constants[92]))))*(states[1]-constants[38]*exp((-2.00000*(algebraic[19]-0.00000))/constants[92])))*states[17])*states[19])*states[18]
    algebraic[18] = constants[92]*log(constants[40]/states[15])
    algebraic[91] = (((constants[86]*constants[73])*(power(1.00000+power(constants[71]/constants[39], 1.20000), -1.00000)))*(power(1.00000+power(constants[72]/states[15], 1.30000), -1.00000)))*(power(1.00000+exp(-((algebraic[19]-algebraic[18])+110.000)/20.0000), -1.00000))
    algebraic[73] = constants[92]*log((constants[40]+0.120000*constants[39])/(states[15]+0.120000*states[14]))
    algebraic[74] = ((constants[56]*(power(states[30], 3.00000)))*states[29])*(algebraic[19]-algebraic[73])
    algebraic[75] = (constants[57]*(power(states[30], 3.00000)))*(algebraic[19]-algebraic[73])
    algebraic[76] = algebraic[74]+algebraic[75]
    algebraic[94] = (states[32]*constants[79])*(algebraic[19]-algebraic[18])
    algebraic[53] = ((((((1.85000e-05*constants[44])*(algebraic[19]-0.00000))/(constants[92]*(1.00000-exp((-1.00000*(algebraic[19]-0.00000))/constants[92]))))*(states[15]-constants[40]*exp((-1.00000*(algebraic[19]-0.00000))/constants[92])))*states[17])*states[19])*states[18]
    algebraic[16] = constants[92]*log(constants[39]/states[14])
    algebraic[69] = ((constants[84]*(algebraic[19]-algebraic[16]))*(0.900000*states[23]+0.100000*states[24]))*states[25]
    algebraic[70] = constants[92]*log((constants[39]+0.120000*constants[40])/(states[14]+0.120000*states[15]))
    algebraic[71] = (constants[85]*(algebraic[19]-algebraic[70]))*(power(states[26], 2.00000))
    algebraic[72] = ((constants[55]*states[27])*states[28])*(algebraic[19]-algebraic[16])
    algebraic[92] = (constants[75]*(algebraic[19]-algebraic[16]))*states[31]
    algebraic[93] = (states[32]*constants[78])*(algebraic[19]-algebraic[16])
    algebraic[36] = ((((((0.000365000*constants[44])*(algebraic[19]-0.00000))/(constants[92]*(1.00000-exp((-1.00000*(algebraic[19]-0.00000))/constants[92]))))*(states[14]-constants[39]*exp((-1.00000*(algebraic[19]-0.00000))/constants[92])))*states[17])*states[19])*states[18]
    algebraic[96] = ((constants[81]*(algebraic[19]-algebraic[16]))*states[33])*states[34]
    algebraic[59] = ((((algebraic[20]+algebraic[36])+algebraic[53])*(1.00000-constants[93]))*1.00000)*constants[94]
    algebraic[67] = custom_piecewise([greater(constants[25] , 0.00000), (((constants[52]*constants[53])*(algebraic[19]-algebraic[16]))*(1.00000+exp((algebraic[19]+20.0000)/20.0000)))*states[22] , True, 0.00000])
    algebraic[95] = algebraic[94]+algebraic[93]
    algebraic[97] = ((((((((((algebraic[95]+algebraic[69])+algebraic[71])+algebraic[96])+algebraic[91])+algebraic[90])+algebraic[76])+algebraic[59])+algebraic[64])+algebraic[67])+algebraic[72])+algebraic[92]
    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)