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 = 121
sizeStates = 31
sizeConstants = 134
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 (second)"
    legend_constants[0] = "protocol in component environment (dimensionless)"
    legend_constants[1] = "L_totmax in component b1_AR_Gs_parameters (uM)"
    legend_constants[2] = "sum_b1_AR in component b1_AR_Gs_parameters (uM)"
    legend_constants[3] = "Gs_tot in component b1_AR_Gs_parameters (uM)"
    legend_constants[4] = "Kl in component b1_AR_Gs_parameters (uM)"
    legend_constants[5] = "Kr in component b1_AR_Gs_parameters (uM)"
    legend_constants[6] = "Kc in component b1_AR_Gs_parameters (uM)"
    legend_constants[7] = "k_bar_kp in component b1_AR_Gs_parameters (per_sec)"
    legend_constants[8] = "k_bar_km in component b1_AR_Gs_parameters (per_sec)"
    legend_constants[9] = "k_p_kap in component b1_AR_Gs_parameters (per_uM_per_sec)"
    legend_constants[10] = "k_p_kam in component b1_AR_Gs_parameters (per_sec)"
    legend_constants[11] = "k_g_act in component b1_AR_Gs_parameters (per_sec)"
    legend_constants[12] = "k_hyd in component b1_AR_Gs_parameters (per_sec)"
    legend_constants[13] = "k_reassoc in component b1_AR_Gs_parameters (per_uM_per_sec)"
    legend_constants[14] = "AC_tot in component cAMP_parameters (uM)"
    legend_constants[15] = "ATP in component cAMP_parameters (uM)"
    legend_constants[16] = "PDE_tot in component cAMP_parameters (uM)"
    legend_constants[17] = "IBMX_tot in component cAMP_parameters (uM)"
    legend_constants[18] = "Fsk_tot in component cAMP_parameters (uM)"
    legend_constants[19] = "k_ac_basal in component cAMP_parameters (per_sec)"
    legend_constants[20] = "k_ac_gsa in component cAMP_parameters (per_sec)"
    legend_constants[21] = "k_ac_fsk in component cAMP_parameters (per_sec)"
    legend_constants[22] = "k_pde in component cAMP_parameters (per_sec)"
    legend_constants[23] = "Km_basal in component cAMP_parameters (uM)"
    legend_constants[24] = "Km_gsa in component cAMP_parameters (uM)"
    legend_constants[25] = "Km_fsk in component cAMP_parameters (uM)"
    legend_constants[26] = "Km_pde in component cAMP_parameters (uM)"
    legend_constants[27] = "K_gsa in component cAMP_parameters (uM)"
    legend_constants[28] = "K_fsk in component cAMP_parameters (uM)"
    legend_constants[29] = "Ki_ibmx in component cAMP_parameters (uM)"
    legend_constants[30] = "PKAI_tot in component PKA_parameters (uM)"
    legend_constants[31] = "PKAII_tot in component PKA_parameters (uM)"
    legend_constants[32] = "PKI_tot in component PKA_parameters (uM)"
    legend_constants[33] = "K_a in component PKA_parameters (uM)"
    legend_constants[34] = "K_b in component PKA_parameters (uM)"
    legend_constants[35] = "K_d in component PKA_parameters (uM)"
    legend_constants[36] = "Ki_pki in component PKA_parameters (uM)"
    legend_constants[37] = "PLB_tot in component PLB_parameters (uM)"
    legend_constants[38] = "PP1_tot in component PLB_parameters (uM)"
    legend_constants[39] = "Inhib1_tot in component PLB_parameters (uM)"
    legend_constants[40] = "k_pka_plb in component PLB_parameters (per_sec)"
    legend_constants[41] = "Km_pka_plb in component PLB_parameters (uM)"
    legend_constants[42] = "k_pp1_plb in component PLB_parameters (per_sec)"
    legend_constants[43] = "Km_pp1_plb in component PLB_parameters (uM)"
    legend_constants[44] = "k_pka_i1 in component PLB_parameters (per_sec)"
    legend_constants[45] = "Km_pka_i1 in component PLB_parameters (uM)"
    legend_constants[46] = "Vmax_pp2a_i1 in component PLB_parameters (uM_per_sec)"
    legend_constants[47] = "Km_pp2a_i1 in component PLB_parameters (uM)"
    legend_constants[48] = "Ki_inhib1 in component PLB_parameters (uM)"
    legend_constants[49] = "epsilon in component LCC_parameters (dimensionless)"
    legend_constants[50] = "LCC_tot in component LCC_parameters (uM)"
    legend_constants[51] = "PP1_lcc_tot in component LCC_parameters (uM)"
    legend_constants[52] = "PP2A_lcc_tot in component LCC_parameters (uM)"
    legend_constants[53] = "k_pka_lcc in component LCC_parameters (per_sec)"
    legend_constants[54] = "Km_pka_lcc in component LCC_parameters (uM)"
    legend_constants[55] = "k_pp1_lcc in component LCC_parameters (per_sec)"
    legend_constants[56] = "Km_pp1_lcc in component LCC_parameters (uM)"
    legend_constants[57] = "k_pp2a_lcc in component LCC_parameters (per_sec)"
    legend_constants[58] = "Km_pp2a_lcc in component LCC_parameters (uM)"
    legend_constants[59] = "V_myo in component EC_Coupling_Parameters (uL)"
    legend_constants[60] = "V_nsr in component EC_Coupling_Parameters (uL)"
    legend_constants[61] = "V_jsr in component EC_Coupling_Parameters (uL)"
    legend_constants[62] = "A_Cap in component EC_Coupling_Parameters (cm2)"
    legend_constants[63] = "Temp in component EC_Coupling_Parameters (kelvin)"
    legend_constants[64] = "Na_ext in component EC_Coupling_Parameters (mM)"
    legend_constants[65] = "K_ext in component EC_Coupling_Parameters (mM)"
    legend_constants[66] = "Ca_ext in component EC_Coupling_Parameters (mM)"
    legend_constants[67] = "G_Na in component EC_Coupling_Parameters (mS_per_uF)"
    legend_constants[68] = "G_to in component EC_Coupling_Parameters (mS_per_uF)"
    legend_constants[69] = "G_ss in component EC_Coupling_Parameters (mS_per_uF)"
    legend_constants[70] = "G_Ki_bar in component EC_Coupling_Parameters (mS_per_uF)"
    legend_constants[71] = "G_Kp in component EC_Coupling_Parameters (mS_per_uF)"
    legend_constants[72] = "f in component EC_Coupling_Parameters (per_sec)"
    legend_constants[73] = "g in component EC_Coupling_Parameters (per_sec)"
    legend_constants[74] = "gamma_o in component EC_Coupling_Parameters (per_mM_per_sec)"
    legend_constants[75] = "omega in component EC_Coupling_Parameters (per_sec)"
    legend_constants[76] = "p_Ca in component EC_Coupling_Parameters (cm_per_sec)"
    legend_constants[77] = "p_K in component EC_Coupling_Parameters (cm_per_sec)"
    legend_constants[78] = "N_lcc in component EC_Coupling_Parameters (dimensionless)"
    legend_constants[79] = "I_Ca05 in component EC_Coupling_Parameters (uA_per_uF)"
    legend_constants[80] = "k_NaCa in component EC_Coupling_Parameters (uA_per_uF)"
    legend_constants[81] = "Km_Na in component EC_Coupling_Parameters (mM)"
    legend_constants[82] = "Km_Ca in component EC_Coupling_Parameters (mM)"
    legend_constants[83] = "k_sat in component EC_Coupling_Parameters (dimensionless)"
    legend_constants[84] = "eta in component EC_Coupling_Parameters (dimensionless)"
    legend_constants[85] = "I_bar_NaK in component EC_Coupling_Parameters (uA_per_uF)"
    legend_constants[86] = "Km_Nai in component EC_Coupling_Parameters (mM)"
    legend_constants[87] = "Km_Ko in component EC_Coupling_Parameters (mM)"
    legend_constants[88] = "I_bar_PCa in component EC_Coupling_Parameters (uA_per_uF)"
    legend_constants[89] = "Km_PCa in component EC_Coupling_Parameters (mM)"
    legend_constants[90] = "G_CaB in component EC_Coupling_Parameters (uA_per_uF)"
    legend_constants[91] = "G_NaB in component EC_Coupling_Parameters (uA_per_uF)"
    legend_constants[92] = "Pns in component EC_Coupling_Parameters (dimensionless)"
    legend_constants[93] = "Km_NS in component EC_Coupling_Parameters (mM)"
    legend_constants[94] = "I_up_bar in component EC_Coupling_Parameters (mM_per_sec)"
    legend_constants[95] = "Km_up0 in component EC_Coupling_Parameters (mM)"
    legend_constants[96] = "NSR_bar in component EC_Coupling_Parameters (mM)"
    legend_constants[97] = "tau_on in component EC_Coupling_Parameters (second)"
    legend_constants[98] = "tau_off in component EC_Coupling_Parameters (second)"
    legend_constants[99] = "G_max_rel in component EC_Coupling_Parameters (mM_per_sec)"
    legend_constants[100] = "d_Cai_th in component EC_Coupling_Parameters (mM)"
    legend_constants[101] = "Km_rel in component EC_Coupling_Parameters (mM)"
    legend_constants[102] = "CSQN_th in component EC_Coupling_Parameters (mM)"
    legend_constants[103] = "CSQN_bar in component EC_Coupling_Parameters (mM)"
    legend_constants[104] = "Km_CSQN in component EC_Coupling_Parameters (mM)"
    legend_constants[105] = "tau_tr in component EC_Coupling_Parameters (second)"
    legend_constants[106] = "TRPN_bar in component EC_Coupling_Parameters (mM)"
    legend_constants[107] = "CMDN_bar in component EC_Coupling_Parameters (mM)"
    legend_constants[108] = "INDO_bar in component EC_Coupling_Parameters (mM)"
    legend_constants[109] = "Km_TRPN in component EC_Coupling_Parameters (mM)"
    legend_constants[110] = "Km_CMDN in component EC_Coupling_Parameters (mM)"
    legend_constants[111] = "Km_INDO in component EC_Coupling_Parameters (mM)"
    legend_algebraic[61] = "LR in component b1_AR_module (uM)"
    legend_algebraic[62] = "LRG in component b1_AR_module (uM)"
    legend_algebraic[63] = "RG in component b1_AR_module (uM)"
    legend_algebraic[79] = "BARK_DESENS in component b1_AR_module (uM_per_sec)"
    legend_algebraic[11] = "BARK_RESENS in component b1_AR_module (uM_per_sec)"
    legend_algebraic[87] = "PKA_DESENS in component b1_AR_module (uM_per_sec)"
    legend_algebraic[22] = "PKA_RESENS in component b1_AR_module (uM_per_sec)"
    legend_algebraic[80] = "G_ACT in component b1_AR_module (uM_per_sec)"
    legend_algebraic[24] = "HYD in component b1_AR_module (uM_per_sec)"
    legend_algebraic[26] = "REASSOC in component b1_AR_module (uM_per_sec)"
    legend_algebraic[64] = "L in component b1_AR_module (uM)"
    legend_algebraic[65] = "R in component b1_AR_module (uM)"
    legend_algebraic[66] = "Gs in component b1_AR_module (uM)"
    legend_states[0] = "b1_AR_tot in component b1_AR_module (uM)"
    legend_states[1] = "b1_AR_d in component b1_AR_module (uM)"
    legend_states[2] = "b1_AR_p in component b1_AR_module (uM)"
    legend_states[3] = "Gs_agtp_tot in component b1_AR_module (uM)"
    legend_states[4] = "Gs_agdp in component b1_AR_module (uM)"
    legend_states[5] = "Gs_bg in component b1_AR_module (uM)"
    legend_algebraic[67] = "PKAC_I in component PKA_module (uM)"
    legend_algebraic[68] = "cAMP in component PKA_module (uM)"
    legend_algebraic[41] = "Gsa_GTP in component cAMP_module (uM)"
    legend_algebraic[49] = "Fsk in component cAMP_module (uM)"
    legend_algebraic[42] = "AC in component cAMP_module (uM)"
    legend_constants[131] = "PDE in component cAMP_module (uM)"
    legend_constants[132] = "IBMX in component cAMP_module (uM)"
    legend_states[6] = "cAMP_tot in component cAMP_module (uM)"
    legend_algebraic[43] = "Gsa_GTP_AC in component cAMP_module (uM)"
    legend_algebraic[50] = "Fsk_AC in component cAMP_module (uM)"
    legend_algebraic[47] = "AC_ACT_GSA in component cAMP_module (uM_per_sec)"
    legend_algebraic[45] = "AC_ACT_BASAL in component cAMP_module (uM_per_sec)"
    legend_algebraic[52] = "AC_ACT_FSK in component cAMP_module (uM_per_sec)"
    legend_algebraic[81] = "PDE_ACT in component cAMP_module (uM_per_sec)"
    legend_constants[133] = "PDE_IBMX in component cAMP_module (uM)"
    legend_algebraic[69] = "PKI in component PKA_module (uM)"
    legend_algebraic[70] = "A2RC_I in component PKA_module (uM)"
    legend_algebraic[71] = "A2R_I in component PKA_module (uM)"
    legend_algebraic[72] = "A2RC_II in component PKA_module (uM)"
    legend_algebraic[73] = "A2R_II in component PKA_module (uM)"
    legend_algebraic[74] = "ARC_I in component PKA_module (uM)"
    legend_algebraic[75] = "ARC_II in component PKA_module (uM)"
    legend_algebraic[76] = "PKA_temp in component PKA_module (uM)"
    legend_algebraic[77] = "PKAC_II in component PKA_module (uM)"
    legend_constants[112] = "zed in component PKA_module (dimensionless)"
    legend_constants[113] = "zed1 in component PKA_module (dimensionless)"
    legend_algebraic[28] = "PLB in component PLB_module (uM)"
    legend_algebraic[82] = "PLB_PHOSPH in component PLB_module (uM_per_sec)"
    legend_algebraic[59] = "PLB_DEPHOSPH in component PLB_module (uM_per_sec)"
    legend_algebraic[29] = "Inhib1 in component PLB_module (uM)"
    legend_algebraic[54] = "Inhib1p_PP1 in component PLB_module (uM)"
    legend_algebraic[83] = "Inhib1_PHOSPH in component PLB_module (uM_per_sec)"
    legend_algebraic[31] = "Inhib1_DEPHOSPH in component PLB_module (uM_per_sec)"
    legend_states[7] = "PLB_p in component PLB_module (uM)"
    legend_states[8] = "Inhib1p_tot in component PLB_module (uM)"
    legend_algebraic[55] = "Inhib1p in component PLB_module (uM)"
    legend_algebraic[56] = "PP1 in component PLB_module (uM)"
    legend_algebraic[0] = "frac_PLB_p in component PLB_module (dimensionless)"
    legend_algebraic[30] = "frac_PLB in component PLB_module (dimensionless)"
    legend_constants[124] = "frac_PLB_o in component PLB_module (dimensionless)"
    legend_algebraic[33] = "LCCa in component LCC_module (uM)"
    legend_algebraic[84] = "LCCa_PHOSPH in component LCC_module (uM_per_sec)"
    legend_algebraic[35] = "LCCa_DEPHOSPH in component LCC_module (uM_per_sec)"
    legend_algebraic[37] = "LCCb in component LCC_module (uM)"
    legend_algebraic[85] = "LCCb_PHOSPH in component LCC_module (uM_per_sec)"
    legend_algebraic[39] = "LCCb_DEPHOSPH in component LCC_module (uM_per_sec)"
    legend_states[9] = "LCCa_p in component LCC_module (uM)"
    legend_states[10] = "LCCb_p in component LCC_module (uM)"
    legend_algebraic[1] = "frac_LCCa_p in component LCC_module (dimensionless)"
    legend_constants[125] = "frac_LCCa_po in component LCC_module (dimensionless)"
    legend_algebraic[32] = "frac_LCCb_p in component LCC_module (dimensionless)"
    legend_constants[126] = "frac_LCCb_po in component LCC_module (dimensionless)"
    legend_algebraic[34] = "E_Na in component Nernst_Potentials (mV)"
    legend_algebraic[36] = "E_K in component Nernst_Potentials (mV)"
    legend_algebraic[38] = "E_Ca in component Nernst_Potentials (mV)"
    legend_constants[127] = "E_Cl in component Nernst_Potentials (mV)"
    legend_constants[114] = "R in component Nernst_Potentials (joules_per_kmole_kelvin)"
    legend_constants[115] = "Frdy in component Nernst_Potentials (coulombs_per_mole)"
    legend_constants[128] = "FoRT in component Nernst_Potentials (per_mV)"
    legend_constants[116] = "z_Na in component Nernst_Potentials (dimensionless)"
    legend_constants[117] = "z_K in component Nernst_Potentials (dimensionless)"
    legend_constants[118] = "z_Ca in component Nernst_Potentials (dimensionless)"
    legend_states[11] = "Na_i in component Ion_Concentrations_and_Membrane_Potential (mM)"
    legend_states[12] = "K_i in component Ion_Concentrations_and_Membrane_Potential (mM)"
    legend_states[13] = "Ca_i in component Ion_Concentrations_and_Membrane_Potential (mM)"
    legend_algebraic[2] = "am in component Fast_Na_Current (per_sec)"
    legend_algebraic[12] = "bm in component Fast_Na_Current (per_sec)"
    legend_algebraic[3] = "ah in component Fast_Na_Current (per_sec)"
    legend_algebraic[4] = "aj in component Fast_Na_Current (per_sec)"
    legend_algebraic[13] = "bh in component Fast_Na_Current (per_sec)"
    legend_algebraic[14] = "bj in component Fast_Na_Current (per_sec)"
    legend_states[14] = "m in component Fast_Na_Current (dimensionless)"
    legend_states[15] = "h in component Fast_Na_Current (dimensionless)"
    legend_states[16] = "j in component Fast_Na_Current (dimensionless)"
    legend_algebraic[40] = "I_Na in component Fast_Na_Current (uA_per_uF)"
    legend_states[17] = "V_m in component Ion_Concentrations_and_Membrane_Potential (mV)"
    legend_algebraic[5] = "a_lcc in component L_Type_Calcium_Current (per_sec)"
    legend_algebraic[15] = "b_lcc in component L_Type_Calcium_Current (per_sec)"
    legend_algebraic[17] = "f_lcc in component L_Type_Calcium_Current (per_sec)"
    legend_algebraic[6] = "y_lcc_inf in component L_Type_Calcium_Current (dimensionless)"
    legend_algebraic[16] = "tau_y_lcc in component L_Type_Calcium_Current (second)"
    legend_algebraic[23] = "gamma in component L_Type_Calcium_Current (per_sec)"
    legend_algebraic[25] = "v_gamma in component L_Type_Calcium_Current (per_sec)"
    legend_algebraic[27] = "v_omega in component L_Type_Calcium_Current (per_sec)"
    legend_states[18] = "v in component L_Type_Calcium_Current (dimensionless)"
    legend_states[19] = "w in component L_Type_Calcium_Current (dimensionless)"
    legend_states[20] = "x in component L_Type_Calcium_Current (dimensionless)"
    legend_states[21] = "y in component L_Type_Calcium_Current (dimensionless)"
    legend_states[22] = "z in component L_Type_Calcium_Current (dimensionless)"
    legend_algebraic[44] = "i_bar_Ca in component L_Type_Calcium_Current (uA_per_uF)"
    legend_algebraic[46] = "i_bar_K in component L_Type_Calcium_Current (uA_per_uF)"
    legend_algebraic[48] = "f_avail in component L_Type_Calcium_Current (dimensionless)"
    legend_algebraic[51] = "I_Ca in component L_Type_Calcium_Current (uA_per_uF)"
    legend_algebraic[53] = "I_CaK in component L_Type_Calcium_Current (uA_per_uF)"
    legend_algebraic[57] = "I_Ca_tot in component L_Type_Calcium_Current (uA_per_uF)"
    legend_constants[119] = "alpha in component L_Type_Calcium_Current (ibar_to_I)"
    legend_algebraic[7] = "r_toss in component Transient_Outward_K_Current (dimensionless)"
    legend_algebraic[8] = "s_toss in component Transient_Outward_K_Current (dimensionless)"
    legend_algebraic[18] = "tau_r_to in component Transient_Outward_K_Current (second)"
    legend_algebraic[19] = "tau_s_to in component Transient_Outward_K_Current (second)"
    legend_algebraic[20] = "tau_ss_to in component Transient_Outward_K_Current (second)"
    legend_states[23] = "r_to in component Transient_Outward_K_Current (dimensionless)"
    legend_states[24] = "s_to in component Transient_Outward_K_Current (dimensionless)"
    legend_states[25] = "ss_to in component Transient_Outward_K_Current (dimensionless)"
    legend_algebraic[58] = "I_to in component Transient_Outward_K_Current (uA_per_uF)"
    legend_algebraic[9] = "r_ss_inf in component Steady_State_K_Current (dimensionless)"
    legend_algebraic[21] = "tau_r_ss in component Steady_State_K_Current (second)"
    legend_algebraic[10] = "s_ss_inf in component Steady_State_K_Current (dimensionless)"
    legend_constants[129] = "tau_s_ss in component Steady_State_K_Current (second)"
    legend_states[26] = "r_ss in component Steady_State_K_Current (dimensionless)"
    legend_states[27] = "s_ss in component Steady_State_K_Current (dimensionless)"
    legend_algebraic[60] = "I_ss in component Steady_State_K_Current (uA_per_uF)"
    legend_algebraic[78] = "a_Ki in component Time_Independent_K_Current (dimensionless)"
    legend_algebraic[86] = "b_Ki in component Time_Independent_K_Current (dimensionless)"
    legend_algebraic[88] = "Ki_ss in component Time_Independent_K_Current (dimensionless)"
    legend_algebraic[89] = "I_Ki in component Time_Independent_K_Current (uA_per_uF)"
    legend_algebraic[90] = "Kp in component Plateau_K_Current (dimensionless)"
    legend_algebraic[91] = "I_Kp in component Plateau_K_Current (uA_per_uF)"
    legend_algebraic[92] = "s4 in component Na_Ca_Exchanger_Current (dimensionless)"
    legend_algebraic[93] = "s5 in component Na_Ca_Exchanger_Current (dimensionless)"
    legend_algebraic[94] = "I_NCX in component Na_Ca_Exchanger_Current (uA_per_uF)"
    legend_constants[130] = "sigma in component Na_K_Pump_Current (dimensionless)"
    legend_algebraic[95] = "f_NaK in component Na_K_Pump_Current (dimensionless)"
    legend_algebraic[96] = "I_NaK in component Na_K_Pump_Current (uA_per_uF)"
    legend_algebraic[97] = "I_PCa in component Sarcolemmal_Ca_Pump_Current (uA_per_uF)"
    legend_algebraic[98] = "I_CaB in component Ca_Background_Current (uA_per_uF)"
    legend_algebraic[99] = "I_NaB in component Na_Background_Current (uA_per_uF)"
    legend_algebraic[100] = "I_Na_tot in component Total_Membrane_Currents (uA_per_uF)"
    legend_algebraic[101] = "I_K_tot in component Total_Membrane_Currents (uA_per_uF)"
    legend_algebraic[102] = "I_Ca_tot in component Total_Membrane_Currents (uA_per_uF)"
    legend_algebraic[103] = "t_rel_sub in component Calcium_Induced_Calcium_Release (second)"
    legend_algebraic[105] = "ryr_on in component Calcium_Induced_Calcium_Release (dimensionless)"
    legend_algebraic[107] = "ryr_off in component Calcium_Induced_Calcium_Release (dimensionless)"
    legend_algebraic[109] = "g_rel in component Calcium_Induced_Calcium_Release (per_sec)"
    legend_algebraic[110] = "I_rel in component Calcium_Induced_Calcium_Release (mM_per_sec)"
    legend_states[28] = "Ca_jsr in component Other_SR_Fluxes_and_Concentrations (mM)"
    legend_states[29] = "trel in component Ion_Concentrations_and_Membrane_Potential (second)"
    legend_algebraic[111] = "Km_up in component Other_SR_Fluxes_and_Concentrations (mM)"
    legend_algebraic[112] = "I_up in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)"
    legend_algebraic[113] = "I_leak in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)"
    legend_algebraic[114] = "I_tr in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)"
    legend_algebraic[116] = "B_jsr in component Other_SR_Fluxes_and_Concentrations (dimensionless)"
    legend_states[30] = "Ca_nsr in component Other_SR_Fluxes_and_Concentrations (mM)"
    legend_algebraic[118] = "SR_content in component Other_SR_Fluxes_and_Concentrations (uM)"
    legend_algebraic[115] = "b_trpn in component Cytoplasmic_Calcium_Buffering (dimensionless)"
    legend_algebraic[117] = "b_cmdn in component Cytoplasmic_Calcium_Buffering (dimensionless)"
    legend_algebraic[119] = "b_indo in component Cytoplasmic_Calcium_Buffering (dimensionless)"
    legend_algebraic[120] = "B_myo in component Cytoplasmic_Calcium_Buffering (dimensionless)"
    legend_algebraic[108] = "I_app in component Ion_Concentrations_and_Membrane_Potential (uA_per_uF)"
    legend_constants[120] = "V_hold in component Ion_Concentrations_and_Membrane_Potential (mV)"
    legend_constants[121] = "V_test in component Ion_Concentrations_and_Membrane_Potential (mV)"
    legend_algebraic[104] = "V_clamp in component Ion_Concentrations_and_Membrane_Potential (mV)"
    legend_constants[122] = "R_clamp in component Ion_Concentrations_and_Membrane_Potential (ms)"
    legend_constants[123] = "lambda in component Ion_Concentrations_and_Membrane_Potential (nFC_per_sAcm2)"
    legend_algebraic[106] = "I_pace in component Ion_Concentrations_and_Membrane_Potential (uA_per_uF)"
    legend_rates[0] = "d/dt b1_AR_tot in component b1_AR_module (uM)"
    legend_rates[1] = "d/dt b1_AR_d in component b1_AR_module (uM)"
    legend_rates[2] = "d/dt b1_AR_p in component b1_AR_module (uM)"
    legend_rates[3] = "d/dt Gs_agtp_tot in component b1_AR_module (uM)"
    legend_rates[4] = "d/dt Gs_agdp in component b1_AR_module (uM)"
    legend_rates[5] = "d/dt Gs_bg in component b1_AR_module (uM)"
    legend_rates[6] = "d/dt cAMP_tot in component cAMP_module (uM)"
    legend_rates[7] = "d/dt PLB_p in component PLB_module (uM)"
    legend_rates[8] = "d/dt Inhib1p_tot in component PLB_module (uM)"
    legend_rates[9] = "d/dt LCCa_p in component LCC_module (uM)"
    legend_rates[10] = "d/dt LCCb_p in component LCC_module (uM)"
    legend_rates[14] = "d/dt m in component Fast_Na_Current (dimensionless)"
    legend_rates[15] = "d/dt h in component Fast_Na_Current (dimensionless)"
    legend_rates[16] = "d/dt j in component Fast_Na_Current (dimensionless)"
    legend_rates[18] = "d/dt v in component L_Type_Calcium_Current (dimensionless)"
    legend_rates[19] = "d/dt w in component L_Type_Calcium_Current (dimensionless)"
    legend_rates[20] = "d/dt x in component L_Type_Calcium_Current (dimensionless)"
    legend_rates[21] = "d/dt y in component L_Type_Calcium_Current (dimensionless)"
    legend_rates[22] = "d/dt z in component L_Type_Calcium_Current (dimensionless)"
    legend_rates[23] = "d/dt r_to in component Transient_Outward_K_Current (dimensionless)"
    legend_rates[24] = "d/dt s_to in component Transient_Outward_K_Current (dimensionless)"
    legend_rates[25] = "d/dt ss_to in component Transient_Outward_K_Current (dimensionless)"
    legend_rates[26] = "d/dt r_ss in component Steady_State_K_Current (dimensionless)"
    legend_rates[27] = "d/dt s_ss in component Steady_State_K_Current (dimensionless)"
    legend_rates[30] = "d/dt Ca_nsr in component Other_SR_Fluxes_and_Concentrations (mM)"
    legend_rates[28] = "d/dt Ca_jsr in component Other_SR_Fluxes_and_Concentrations (mM)"
    legend_rates[11] = "d/dt Na_i in component Ion_Concentrations_and_Membrane_Potential (mM)"
    legend_rates[12] = "d/dt K_i in component Ion_Concentrations_and_Membrane_Potential (mM)"
    legend_rates[13] = "d/dt Ca_i in component Ion_Concentrations_and_Membrane_Potential (mM)"
    legend_rates[17] = "d/dt V_m in component Ion_Concentrations_and_Membrane_Potential (mV)"
    legend_rates[29] = "d/dt trel in component Ion_Concentrations_and_Membrane_Potential (second)"
    return (legend_states, legend_algebraic, legend_voi, legend_constants)

def initConsts():
    constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
    constants[0] = 1
    constants[1] = 0.01
    constants[2] = 0.0132
    constants[3] = 3.83
    constants[4] = 0.285
    constants[5] = 0.062
    constants[6] = 33
    constants[7] = 1.1e-3
    constants[8] = 2.2e-3
    constants[9] = 3.6e-3
    constants[10] = 2.2e-3
    constants[11] = 16
    constants[12] = 0.8
    constants[13] = 1.21e3
    constants[14] = 49.7e-3
    constants[15] = 5e3
    constants[16] = 38.9e-3
    constants[17] = 0
    constants[18] = 0
    constants[19] = 0.2
    constants[20] = 8.5
    constants[21] = 7.3
    constants[22] = 5
    constants[23] = 1.03e3
    constants[24] = 315
    constants[25] = 860
    constants[26] = 1.3
    constants[27] = 0.4
    constants[28] = 44
    constants[29] = 30
    constants[30] = 0.59
    constants[31] = 0.025
    constants[32] = 0.18
    constants[33] = 9.14
    constants[34] = 1.64
    constants[35] = 4.375
    constants[36] = 0.2e-3
    constants[37] = 106
    constants[38] = 0.89
    constants[39] = 0.3
    constants[40] = 54
    constants[41] = 21
    constants[42] = 8.5
    constants[43] = 7
    constants[44] = 60
    constants[45] = 1
    constants[46] = 14
    constants[47] = 1
    constants[48] = 1e-3
    constants[49] = 10
    constants[50] = 25e-3
    constants[51] = 25e-3
    constants[52] = 25e-3
    constants[53] = 54
    constants[54] = 21
    constants[55] = 8.52
    constants[56] = 3
    constants[57] = 10.1
    constants[58] = 3
    constants[59] = 20.8e-6
    constants[60] = 9.88e-7
    constants[61] = 9.3e-8
    constants[62] = 1.534e-4
    constants[63] = 310
    constants[64] = 140
    constants[65] = 5.4
    constants[66] = 1.8
    constants[67] = 8
    constants[68] = 0.35
    constants[69] = 0.07
    constants[70] = 0.24
    constants[71] = 0.008
    constants[72] = 300
    constants[73] = 2e3
    constants[74] = 5187.5
    constants[75] = 10
    constants[76] = 1.7469e-8
    constants[77] = 3.234e-11
    constants[78] = 3e5
    constants[79] = -0.458
    constants[80] = 1483
    constants[81] = 87.5
    constants[82] = 1.38
    constants[83] = 0.1
    constants[84] = 0.35
    constants[85] = 1.1
    constants[86] = 10
    constants[87] = 1.5
    constants[88] = 1.15
    constants[89] = 0.5e-3
    constants[90] = 2.8e-3
    constants[91] = 1.18e-3
    constants[92] = 0
    constants[93] = 1.2e-3
    constants[94] = 4.7
    constants[95] = 3e-4
    constants[96] = 15
    constants[97] = 2e-3
    constants[98] = 2e-3
    constants[99] = 60e3
    constants[100] = 0.18e-3
    constants[101] = 0.8e-3
    constants[102] = 8.75
    constants[103] = 15
    constants[104] = 0.8
    constants[105] = 5.7e-4
    constants[106] = 0.07
    constants[107] = 0.05
    constants[108] = 0.07
    constants[109] = 0.5128e-3
    constants[110] = 2.38e-3
    constants[111] = 8.44e-4
    states[0] = 0.01205
    states[1] = 0
    states[2] = 1.154e-3
    states[3] = 0.02505
    states[4] = 6.446e-4
    states[5] = 0.02569
    states[6] = 0.8453
    constants[112] = 0
    constants[113] = 0
    states[7] = 4.105
    states[8] = 0.0526
    states[9] = 5.103e-3
    states[10] = 5.841e-3
    constants[114] = 8314
    constants[115] = 96485
    constants[116] = 1
    constants[117] = 1
    constants[118] = 2
    states[11] = 16
    states[12] = 145
    states[13] = 1.58e-4
    states[14] = 1.4e-3
    states[15] = 0.99
    states[16] = 0.99
    states[17] = -85.66
    states[18] = 0
    states[19] = 0
    states[20] = 0.13
    states[21] = 0.96
    states[22] = 0.92
    constants[119] = 10e5
    states[23] = 1.4e-3
    states[24] = 1
    states[25] = 0.613
    states[26] = 198e-3
    states[27] = 0.43
    states[28] = 1.92
    states[29] = 0.9
    states[30] = 1.92
    constants[120] = -40
    constants[121] = -10
    constants[122] = 0.02
    constants[123] = -1e3
    constants[124] = 0.961300
    constants[125] = 0.204100
    constants[126] = 0.233600
    constants[127] = -40.0000
    constants[128] = (constants[115]/constants[114])/constants[63]
    constants[129] = 2.10000
    constants[130] = (exp(constants[64]/67.3000)-1.00000)/7.00000
    rootfind_0(voi, constants, rates, states, algebraic)
    return (states, constants)

def computeRates(voi, states, constants):
    rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
    algebraic[10] = 1.00000/(1.00000+exp((states[17]+87.5000)/10.3000))
    rates[27] = (algebraic[10]-states[27])/constants[129]
    algebraic[2] = (0.320000*(states[17]+47.1300))/(1.00000-exp(-0.100000*(states[17]+47.1300)))
    algebraic[12] = 0.0800000*exp(-states[17]/11.0000)
    rates[14] = 1000.00*(algebraic[2]*(1.00000-states[14])-algebraic[12]*states[14])
    algebraic[3] = custom_piecewise([greater_equal(states[17] , -40.0000), 0.00000 , True, 0.135000*exp((80.0000+states[17])/-6.80000)])
    algebraic[13] = custom_piecewise([greater_equal(states[17] , -40.0000), 1.00000/(0.130000*(1.00000+exp(-(states[17]+10.6600)/11.1000))) , True, 3.56000*exp(0.0790000*states[17])+310000.*exp(0.350000*states[17])])
    rates[15] = 1000.00*(algebraic[3]*(1.00000-states[15])-algebraic[13]*states[15])
    algebraic[4] = custom_piecewise([greater_equal(states[17] , -40.0000), 0.00000 , True, ((-127140.*exp(0.244400*states[17])-3.47400e-05*exp(-0.0439100*states[17]))*(states[17]+37.7800))/(1.00000+exp(0.311000*(states[17]+79.2300)))])
    algebraic[14] = custom_piecewise([greater_equal(states[17] , -40.0000), (0.300000*exp(-2.57500e-07*states[17]))/(1.00000+exp(-0.100000*(states[17]+32.0000))) , True, (0.121200*exp(-0.0105200*states[17]))/(1.00000+exp(-0.137800*(states[17]+40.1400)))])
    rates[16] = 1000.00*(algebraic[4]*(1.00000-states[16])-algebraic[14]*states[16])
    algebraic[5] = 400.000*exp((states[17]+2.00000)/10.0000)
    algebraic[15] = 50.0000*exp((-1.00000*(states[17]+2.00000))/13.0000)
    rates[18] = algebraic[5]*(1.00000-states[18])-algebraic[15]*states[18]
    rates[19] = 2.00000*algebraic[5]*(1.00000-states[19])-(algebraic[15]/2.00000)*states[19]
    algebraic[1] = states[9]/constants[50]
    algebraic[17] = constants[72]*((0.375000*algebraic[1])/constants[125]+0.625000)
    rates[20] = algebraic[17]*(1.00000-states[20])-constants[73]*states[20]
    algebraic[6] = 1.00000/(1.00000+exp((states[17]+55.0000)/7.50000))+0.100000/(1.00000+exp((-states[17]+21.0000)/6.00000))
    algebraic[16] = 0.0200000+0.300000/(1.00000+exp((states[17]+30.0000)/9.50000))
    rates[21] = (algebraic[6]-states[21])/algebraic[16]
    algebraic[7] = 1.00000/(1.00000+exp((states[17]+10.6000)/-11.4200))
    algebraic[18] = 1.00000/(45.1600*exp(0.0357700*(states[17]+50.0000))+98.9000*exp(-0.100000*(states[17]+38.0000)))
    rates[23] = (algebraic[7]-states[23])/algebraic[18]
    algebraic[8] = 1.00000/(1.00000+exp((states[17]+43.5000)/6.88410))
    algebraic[19] = 0.350000*exp(-1.00000*(power((states[17]+70.0000)/15.0000, 2.00000)))+0.0350000
    rates[24] = (algebraic[8]-states[24])/algebraic[19]
    algebraic[20] = 3.70000*exp(-1.00000*(power((states[17]+70.0000)/30.0000, 2.00000)))+0.0350000
    rates[25] = (algebraic[8]-states[25])/algebraic[20]
    algebraic[9] = 1.00000/(1.00000+exp((states[17]+11.5000)/-11.8200))
    algebraic[21] = 10.0000/(45.1600*exp(0.0357700*(states[17]+50.0000))+98.9000*exp(-0.100000*(states[17]+38.0000)))
    rates[26] = (algebraic[9]-states[26])/algebraic[21]
    algebraic[24] = constants[12]*states[3]
    algebraic[26] = constants[13]*states[4]*states[5]
    rates[4] = algebraic[24]-algebraic[26]
    algebraic[23] = constants[74]*states[13]
    algebraic[25] = algebraic[23]*(power(1.00000-states[18], 4.00000)+2.00000*states[18]*(power(1.00000-states[18], 3.00000))+4.00000*(power(states[18], 2.00000))*(power(1.00000-states[18], 2.00000))+8.00000*(power(states[18], 3.00000))*(1.00000-states[18])+16.0000*(power(states[18], 4.00000))*(1.00000-algebraic[17]/constants[73]))
    algebraic[27] = constants[75]*(power(1.00000-states[19], 4.00000)+0.500000*states[19]*(power(1.00000-states[19], 3.00000))+0.250000*(power(states[19], 2.00000))*(power(1.00000-states[19], 2.00000))+0.125000*(power(states[19], 3.00000))*(1.00000-states[19])+0.0625000*(power(states[19], 4.00000)))
    rates[22] = algebraic[27]*(1.00000-states[22])-algebraic[25]*states[22]
    rootfind_1(voi, constants, rates, states, algebraic)
    algebraic[79] = constants[7]*(algebraic[61]+algebraic[62])
    algebraic[11] = constants[8]*states[1]
    rates[1] = algebraic[79]-algebraic[11]
    algebraic[80] = constants[11]*(algebraic[63]+algebraic[62])
    rates[3] = algebraic[80]-algebraic[24]
    rates[5] = algebraic[80]-algebraic[26]
    rootfind_2(voi, constants, rates, states, algebraic)
    algebraic[47] = (constants[20]*algebraic[43]*constants[15])/(constants[24]+constants[15])
    algebraic[45] = (constants[19]*algebraic[42]*constants[15])/(constants[23]+constants[15])
    rootfind_3(voi, constants, rates, states, algebraic)
    algebraic[52] = (constants[21]*algebraic[50]*constants[15])/(constants[25]+constants[15])
    algebraic[81] = (constants[22]*constants[131]*algebraic[68])/(constants[26]+algebraic[68])
    rates[6] = (algebraic[45]+algebraic[47]+algebraic[52])-algebraic[81]
    algebraic[28] = constants[37]-states[7]
    algebraic[82] = (constants[40]*algebraic[67]*algebraic[28])/(constants[41]+algebraic[28])
    rootfind_4(voi, constants, rates, states, algebraic)
    algebraic[59] = (constants[42]*algebraic[56]*states[7])/(constants[43]+states[7])
    rates[7] = algebraic[82]-algebraic[59]
    algebraic[29] = constants[39]-states[8]
    algebraic[83] = (constants[44]*algebraic[67]*algebraic[29])/(constants[45]+algebraic[29])
    algebraic[31] = (constants[46]*states[8])/(constants[47]+states[8])
    rates[8] = algebraic[83]-algebraic[31]
    algebraic[33] = constants[50]-states[9]
    algebraic[84] = (constants[49]*constants[53]*algebraic[77]*algebraic[33])/(constants[54]+constants[49]*algebraic[33])
    algebraic[35] = (constants[49]*constants[57]*constants[52]*states[9])/(constants[58]+constants[49]*states[9])
    rates[9] = algebraic[84]-algebraic[35]
    algebraic[37] = constants[50]-states[10]
    algebraic[85] = (constants[49]*constants[53]*algebraic[77]*algebraic[37])/(constants[54]+constants[49]*algebraic[37])
    algebraic[39] = (constants[49]*constants[55]*constants[51]*states[10])/(constants[56]+constants[49]*states[10])
    rates[10] = algebraic[85]-algebraic[39]
    algebraic[87] = constants[9]*algebraic[67]*states[0]
    algebraic[22] = constants[10]*states[2]
    rates[0] = ((algebraic[11]-algebraic[79])+algebraic[22])-algebraic[87]
    rates[2] = algebraic[87]-algebraic[22]
    algebraic[34] = (1.00000/(constants[128]*constants[116]))*log(constants[64]/states[11])
    algebraic[40] = constants[67]*(power(states[14], 3.00000))*states[15]*states[16]*(states[17]-algebraic[34])
    algebraic[92] = 1.00000*exp(constants[84]*states[17]*constants[128])*(power(states[11], 3.00000))*constants[66]
    algebraic[93] = 1.00000*exp((constants[84]-1.00000)*states[17]*constants[128])*(power(constants[64], 3.00000))*states[13]
    algebraic[94] = (constants[80]/((power(constants[81], 3.00000)+power(constants[64], 3.00000))*(constants[82]+constants[66])*(1.00000+constants[83]*exp((constants[84]-1.00000)*states[17]*constants[128]))))*(algebraic[92]-algebraic[93])
    algebraic[95] = 1.00000/(1.00000+0.124500*exp(-0.100000*states[17]*constants[128])+0.0365000*constants[130]*exp(-states[17]*constants[128]))
    algebraic[96] = (((constants[85]*algebraic[95])/(1.00000+power(constants[86]/states[11], 1.50000)))*constants[65])/(constants[65]+constants[87])
    algebraic[99] = constants[91]*(states[17]-algebraic[34])
    algebraic[100] = algebraic[40]+algebraic[99]+3.00000*algebraic[94]+3.00000*algebraic[96]
    rates[11] = (constants[123]*algebraic[100]*constants[62])/(constants[59]*constants[116]*constants[115])
    algebraic[46] = (constants[119]*constants[77]*states[17]*constants[115]*constants[128]*(states[12]*exp(states[17]*constants[128])-constants[65]))/(exp(states[17]*constants[128])-1.00000)
    algebraic[32] = states[10]/constants[50]
    algebraic[48] = 0.500000*((0.400000*algebraic[32])/constants[126]+0.600000)
    algebraic[44] = (constants[119]*constants[76]*4.00000*states[17]*constants[115]*constants[128]*(0.00100000*exp(2.00000*states[17]*constants[128])-0.341000*constants[66]))/(exp(2.00000*states[17]*constants[128])-1.00000)
    algebraic[51] = algebraic[44]*constants[78]*algebraic[48]*(power(states[18], 4.00000))*states[20]*states[21]*states[22]
    algebraic[53] = (algebraic[46]/(1.00000+algebraic[51]/constants[79]))*constants[78]*algebraic[48]*(power(states[18], 4.00000))*states[20]*states[21]*states[22]
    algebraic[36] = (1.00000/(constants[128]*constants[117]))*log(constants[65]/states[12])
    algebraic[58] = constants[68]*states[23]*(0.886000*states[24]+0.114000*states[25])*(states[17]-algebraic[36])
    algebraic[60] = constants[69]*states[26]*states[27]*(states[17]-algebraic[36])
    algebraic[78] = 1.02000/(1.00000+exp(0.238500*((states[17]-algebraic[36])-59.2150)))
    algebraic[86] = (0.491240*exp(0.0803200*((states[17]+5.47600)-algebraic[36]))+exp(0.0617500*((states[17]-algebraic[36])-594.310)))/(1.00000+exp(-0.514300*((states[17]-algebraic[36])+4.75300)))
    algebraic[88] = algebraic[78]/(algebraic[78]+algebraic[86])
    algebraic[89] = constants[70]*(power(constants[65]/5.40000, 1.0/2))*algebraic[88]*(states[17]-algebraic[36])
    algebraic[90] = 1.00000/(1.00000+exp((7.48800-states[17])/5.98000))
    algebraic[91] = constants[71]*algebraic[90]*(states[17]-algebraic[36])
    algebraic[101] = ((algebraic[58]+algebraic[60]+algebraic[89]+algebraic[91])-2.00000*algebraic[96])+algebraic[53]
    rates[12] = (constants[123]*algebraic[101]*constants[62])/(constants[59]*constants[117]*constants[115])
    algebraic[97] = (constants[88]*states[13])/(constants[89]+states[13])
    algebraic[38] = (1.00000/(constants[128]*constants[118]))*log(constants[66]/states[13])
    algebraic[98] = constants[90]*(states[17]-algebraic[38])
    algebraic[102] = (algebraic[51]+algebraic[98]+algebraic[97])-2.00000*algebraic[94]
    algebraic[104] = custom_piecewise([greater(voi , 59.1000) & less(voi , 59.5000), constants[121] , True, constants[120]])
    algebraic[106] = custom_piecewise([greater(sin(2.00000* pi*voi) , 0.999900), 10.0000 , True, 0.00000])
    algebraic[108] = custom_piecewise([equal(constants[0] , 0.00000), 0.00000 , equal(constants[0] , 1.00000), algebraic[106] , True, (algebraic[104]-states[17])/constants[122]])
    rates[17] = -1000.00*((algebraic[102]+algebraic[101]+algebraic[100])-algebraic[108])
    rates[29] = custom_piecewise([greater(-1000.00*((algebraic[102]+algebraic[101]+algebraic[100])-algebraic[108]) , 30000.0), 1.00000-10000.0*states[29] , True, 1.00000])
    algebraic[30] = algebraic[28]/constants[37]
    algebraic[111] = (constants[95]*(1.00000+2.00000*algebraic[30]))/(1.00000+2.00000*constants[124])
    algebraic[112] = (constants[94]*(power(states[13], 2.00000)))/(power(algebraic[111], 2.00000)+power(states[13], 2.00000))
    algebraic[113] = (constants[94]*states[30])/constants[96]
    algebraic[114] = (states[30]-states[28])/constants[105]
    rates[30] = (algebraic[112]-algebraic[113])-(algebraic[114]*constants[61])/constants[60]
    algebraic[103] = states[29]+0.00200000
    algebraic[105] = 1.00000-exp(-algebraic[103]/constants[97])
    algebraic[107] = exp(-algebraic[103]/constants[98])
    algebraic[109] = constants[99]/(1.00000*(1.00000+exp((algebraic[102]+5.00000)/0.900000)))
    algebraic[110] = algebraic[109]*algebraic[105]*algebraic[107]*(states[28]-states[13])
    algebraic[116] = 1.00000/(1.00000+(constants[103]*constants[104])/(power(constants[104]+states[28], 2.00000)))
    rates[28] = algebraic[116]*(algebraic[114]-algebraic[110])
    algebraic[115] = (constants[106]*constants[109])/(power(constants[109]+states[13], 2.00000))
    algebraic[117] = (constants[107]*constants[110])/(power(constants[110]+states[13], 2.00000))
    algebraic[119] = (constants[108]*constants[111])/(power(constants[111]+states[13], 2.00000))
    algebraic[120] = 1.00000/(1.00000+algebraic[117]+2.00000*algebraic[115]+algebraic[119])
    rates[13] = algebraic[120]*(((constants[123]*algebraic[102]*constants[62])/(constants[59]*constants[118]*constants[115])-((algebraic[112]-algebraic[113])*constants[60])/constants[59])+(algebraic[110]*constants[61])/constants[59])
    return(rates)

def computeAlgebraic(constants, states, voi):
    algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
    states = array(states)
    voi = array(voi)
    algebraic[10] = 1.00000/(1.00000+exp((states[17]+87.5000)/10.3000))
    algebraic[2] = (0.320000*(states[17]+47.1300))/(1.00000-exp(-0.100000*(states[17]+47.1300)))
    algebraic[12] = 0.0800000*exp(-states[17]/11.0000)
    algebraic[3] = custom_piecewise([greater_equal(states[17] , -40.0000), 0.00000 , True, 0.135000*exp((80.0000+states[17])/-6.80000)])
    algebraic[13] = custom_piecewise([greater_equal(states[17] , -40.0000), 1.00000/(0.130000*(1.00000+exp(-(states[17]+10.6600)/11.1000))) , True, 3.56000*exp(0.0790000*states[17])+310000.*exp(0.350000*states[17])])
    algebraic[4] = custom_piecewise([greater_equal(states[17] , -40.0000), 0.00000 , True, ((-127140.*exp(0.244400*states[17])-3.47400e-05*exp(-0.0439100*states[17]))*(states[17]+37.7800))/(1.00000+exp(0.311000*(states[17]+79.2300)))])
    algebraic[14] = custom_piecewise([greater_equal(states[17] , -40.0000), (0.300000*exp(-2.57500e-07*states[17]))/(1.00000+exp(-0.100000*(states[17]+32.0000))) , True, (0.121200*exp(-0.0105200*states[17]))/(1.00000+exp(-0.137800*(states[17]+40.1400)))])
    algebraic[5] = 400.000*exp((states[17]+2.00000)/10.0000)
    algebraic[15] = 50.0000*exp((-1.00000*(states[17]+2.00000))/13.0000)
    algebraic[1] = states[9]/constants[50]
    algebraic[17] = constants[72]*((0.375000*algebraic[1])/constants[125]+0.625000)
    algebraic[6] = 1.00000/(1.00000+exp((states[17]+55.0000)/7.50000))+0.100000/(1.00000+exp((-states[17]+21.0000)/6.00000))
    algebraic[16] = 0.0200000+0.300000/(1.00000+exp((states[17]+30.0000)/9.50000))
    algebraic[7] = 1.00000/(1.00000+exp((states[17]+10.6000)/-11.4200))
    algebraic[18] = 1.00000/(45.1600*exp(0.0357700*(states[17]+50.0000))+98.9000*exp(-0.100000*(states[17]+38.0000)))
    algebraic[8] = 1.00000/(1.00000+exp((states[17]+43.5000)/6.88410))
    algebraic[19] = 0.350000*exp(-1.00000*(power((states[17]+70.0000)/15.0000, 2.00000)))+0.0350000
    algebraic[20] = 3.70000*exp(-1.00000*(power((states[17]+70.0000)/30.0000, 2.00000)))+0.0350000
    algebraic[9] = 1.00000/(1.00000+exp((states[17]+11.5000)/-11.8200))
    algebraic[21] = 10.0000/(45.1600*exp(0.0357700*(states[17]+50.0000))+98.9000*exp(-0.100000*(states[17]+38.0000)))
    algebraic[24] = constants[12]*states[3]
    algebraic[26] = constants[13]*states[4]*states[5]
    algebraic[23] = constants[74]*states[13]
    algebraic[25] = algebraic[23]*(power(1.00000-states[18], 4.00000)+2.00000*states[18]*(power(1.00000-states[18], 3.00000))+4.00000*(power(states[18], 2.00000))*(power(1.00000-states[18], 2.00000))+8.00000*(power(states[18], 3.00000))*(1.00000-states[18])+16.0000*(power(states[18], 4.00000))*(1.00000-algebraic[17]/constants[73]))
    algebraic[27] = constants[75]*(power(1.00000-states[19], 4.00000)+0.500000*states[19]*(power(1.00000-states[19], 3.00000))+0.250000*(power(states[19], 2.00000))*(power(1.00000-states[19], 2.00000))+0.125000*(power(states[19], 3.00000))*(1.00000-states[19])+0.0625000*(power(states[19], 4.00000)))
    algebraic[79] = constants[7]*(algebraic[61]+algebraic[62])
    algebraic[11] = constants[8]*states[1]
    algebraic[80] = constants[11]*(algebraic[63]+algebraic[62])
    algebraic[47] = (constants[20]*algebraic[43]*constants[15])/(constants[24]+constants[15])
    algebraic[45] = (constants[19]*algebraic[42]*constants[15])/(constants[23]+constants[15])
    algebraic[52] = (constants[21]*algebraic[50]*constants[15])/(constants[25]+constants[15])
    algebraic[81] = (constants[22]*constants[131]*algebraic[68])/(constants[26]+algebraic[68])
    algebraic[28] = constants[37]-states[7]
    algebraic[82] = (constants[40]*algebraic[67]*algebraic[28])/(constants[41]+algebraic[28])
    algebraic[59] = (constants[42]*algebraic[56]*states[7])/(constants[43]+states[7])
    algebraic[29] = constants[39]-states[8]
    algebraic[83] = (constants[44]*algebraic[67]*algebraic[29])/(constants[45]+algebraic[29])
    algebraic[31] = (constants[46]*states[8])/(constants[47]+states[8])
    algebraic[33] = constants[50]-states[9]
    algebraic[84] = (constants[49]*constants[53]*algebraic[77]*algebraic[33])/(constants[54]+constants[49]*algebraic[33])
    algebraic[35] = (constants[49]*constants[57]*constants[52]*states[9])/(constants[58]+constants[49]*states[9])
    algebraic[37] = constants[50]-states[10]
    algebraic[85] = (constants[49]*constants[53]*algebraic[77]*algebraic[37])/(constants[54]+constants[49]*algebraic[37])
    algebraic[39] = (constants[49]*constants[55]*constants[51]*states[10])/(constants[56]+constants[49]*states[10])
    algebraic[87] = constants[9]*algebraic[67]*states[0]
    algebraic[22] = constants[10]*states[2]
    algebraic[34] = (1.00000/(constants[128]*constants[116]))*log(constants[64]/states[11])
    algebraic[40] = constants[67]*(power(states[14], 3.00000))*states[15]*states[16]*(states[17]-algebraic[34])
    algebraic[92] = 1.00000*exp(constants[84]*states[17]*constants[128])*(power(states[11], 3.00000))*constants[66]
    algebraic[93] = 1.00000*exp((constants[84]-1.00000)*states[17]*constants[128])*(power(constants[64], 3.00000))*states[13]
    algebraic[94] = (constants[80]/((power(constants[81], 3.00000)+power(constants[64], 3.00000))*(constants[82]+constants[66])*(1.00000+constants[83]*exp((constants[84]-1.00000)*states[17]*constants[128]))))*(algebraic[92]-algebraic[93])
    algebraic[95] = 1.00000/(1.00000+0.124500*exp(-0.100000*states[17]*constants[128])+0.0365000*constants[130]*exp(-states[17]*constants[128]))
    algebraic[96] = (((constants[85]*algebraic[95])/(1.00000+power(constants[86]/states[11], 1.50000)))*constants[65])/(constants[65]+constants[87])
    algebraic[99] = constants[91]*(states[17]-algebraic[34])
    algebraic[100] = algebraic[40]+algebraic[99]+3.00000*algebraic[94]+3.00000*algebraic[96]
    algebraic[46] = (constants[119]*constants[77]*states[17]*constants[115]*constants[128]*(states[12]*exp(states[17]*constants[128])-constants[65]))/(exp(states[17]*constants[128])-1.00000)
    algebraic[32] = states[10]/constants[50]
    algebraic[48] = 0.500000*((0.400000*algebraic[32])/constants[126]+0.600000)
    algebraic[44] = (constants[119]*constants[76]*4.00000*states[17]*constants[115]*constants[128]*(0.00100000*exp(2.00000*states[17]*constants[128])-0.341000*constants[66]))/(exp(2.00000*states[17]*constants[128])-1.00000)
    algebraic[51] = algebraic[44]*constants[78]*algebraic[48]*(power(states[18], 4.00000))*states[20]*states[21]*states[22]
    algebraic[53] = (algebraic[46]/(1.00000+algebraic[51]/constants[79]))*constants[78]*algebraic[48]*(power(states[18], 4.00000))*states[20]*states[21]*states[22]
    algebraic[36] = (1.00000/(constants[128]*constants[117]))*log(constants[65]/states[12])
    algebraic[58] = constants[68]*states[23]*(0.886000*states[24]+0.114000*states[25])*(states[17]-algebraic[36])
    algebraic[60] = constants[69]*states[26]*states[27]*(states[17]-algebraic[36])
    algebraic[78] = 1.02000/(1.00000+exp(0.238500*((states[17]-algebraic[36])-59.2150)))
    algebraic[86] = (0.491240*exp(0.0803200*((states[17]+5.47600)-algebraic[36]))+exp(0.0617500*((states[17]-algebraic[36])-594.310)))/(1.00000+exp(-0.514300*((states[17]-algebraic[36])+4.75300)))
    algebraic[88] = algebraic[78]/(algebraic[78]+algebraic[86])
    algebraic[89] = constants[70]*(power(constants[65]/5.40000, 1.0/2))*algebraic[88]*(states[17]-algebraic[36])
    algebraic[90] = 1.00000/(1.00000+exp((7.48800-states[17])/5.98000))
    algebraic[91] = constants[71]*algebraic[90]*(states[17]-algebraic[36])
    algebraic[101] = ((algebraic[58]+algebraic[60]+algebraic[89]+algebraic[91])-2.00000*algebraic[96])+algebraic[53]
    algebraic[97] = (constants[88]*states[13])/(constants[89]+states[13])
    algebraic[38] = (1.00000/(constants[128]*constants[118]))*log(constants[66]/states[13])
    algebraic[98] = constants[90]*(states[17]-algebraic[38])
    algebraic[102] = (algebraic[51]+algebraic[98]+algebraic[97])-2.00000*algebraic[94]
    algebraic[104] = custom_piecewise([greater(voi , 59.1000) & less(voi , 59.5000), constants[121] , True, constants[120]])
    algebraic[106] = custom_piecewise([greater(sin(2.00000* pi*voi) , 0.999900), 10.0000 , True, 0.00000])
    algebraic[108] = custom_piecewise([equal(constants[0] , 0.00000), 0.00000 , equal(constants[0] , 1.00000), algebraic[106] , True, (algebraic[104]-states[17])/constants[122]])
    algebraic[30] = algebraic[28]/constants[37]
    algebraic[111] = (constants[95]*(1.00000+2.00000*algebraic[30]))/(1.00000+2.00000*constants[124])
    algebraic[112] = (constants[94]*(power(states[13], 2.00000)))/(power(algebraic[111], 2.00000)+power(states[13], 2.00000))
    algebraic[113] = (constants[94]*states[30])/constants[96]
    algebraic[114] = (states[30]-states[28])/constants[105]
    algebraic[103] = states[29]+0.00200000
    algebraic[105] = 1.00000-exp(-algebraic[103]/constants[97])
    algebraic[107] = exp(-algebraic[103]/constants[98])
    algebraic[109] = constants[99]/(1.00000*(1.00000+exp((algebraic[102]+5.00000)/0.900000)))
    algebraic[110] = algebraic[109]*algebraic[105]*algebraic[107]*(states[28]-states[13])
    algebraic[116] = 1.00000/(1.00000+(constants[103]*constants[104])/(power(constants[104]+states[28], 2.00000)))
    algebraic[115] = (constants[106]*constants[109])/(power(constants[109]+states[13], 2.00000))
    algebraic[117] = (constants[107]*constants[110])/(power(constants[110]+states[13], 2.00000))
    algebraic[119] = (constants[108]*constants[111])/(power(constants[111]+states[13], 2.00000))
    algebraic[120] = 1.00000/(1.00000+algebraic[117]+2.00000*algebraic[115]+algebraic[119])
    algebraic[0] = states[7]/constants[37]
    algebraic[57] = algebraic[51]+algebraic[53]
    algebraic[118] = 1000.00*(((states[28]+states[28]/algebraic[116])*constants[61])/constants[59]+(states[30]*constants[60])/constants[59])
    return algebraic

initialGuess0 = None
def rootfind_0(voi, constants, rates, states, algebraic):
    """Calculate values of algebraic variables for DAE"""
    from scipy.optimize import fsolve
    global initialGuess0
    if initialGuess0 is None: initialGuess0 = ones(3)*0.1
    if not iterable(voi):
        soln = fsolve(residualSN_0, initialGuess0, args=(algebraic, voi, constants, rates, states), xtol=1E-6)
        initialGuess0 = soln
        constants[131] = soln[0]
        constants[132] = soln[1]
        constants[133] = soln[2]
    else:
        for (i,t) in enumerate(voi):
            soln = fsolve(residualSN_0, initialGuess0, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6)
            initialGuess0 = soln
            constants[131][i] = soln[0]
            constants[132][i] = soln[1]
            constants[133][i] = soln[2]

def residualSN_0(algebraicCandidate, algebraic, voi, constants, rates, states):
    resid = array([0.0] * 3)
    constants[131] = algebraicCandidate[0]
    constants[132] = algebraicCandidate[1]
    constants[133] = algebraicCandidate[2]
    resid[0] = (constants[133]-(constants[131]*constants[132])/constants[29])
    resid[1] = (constants[131]-(constants[16]-constants[133]))
    resid[2] = (constants[132]-(constants[17]-constants[133]))
    return resid

initialGuess1 = None
def rootfind_1(voi, constants, rates, states, algebraic):
    """Calculate values of algebraic variables for DAE"""
    from scipy.optimize import fsolve
    global initialGuess1
    if initialGuess1 is None: initialGuess1 = ones(17)*0.1
    if not iterable(voi):
        soln = fsolve(residualSN_1, initialGuess1, args=(algebraic, voi, constants, rates, states), xtol=1E-6)
        initialGuess1 = soln
        algebraic[61] = soln[0]
        algebraic[62] = soln[1]
        algebraic[63] = soln[2]
        algebraic[64] = soln[3]
        algebraic[65] = soln[4]
        algebraic[66] = soln[5]
        algebraic[67] = soln[6]
        algebraic[68] = soln[7]
        algebraic[69] = soln[8]
        algebraic[70] = soln[9]
        algebraic[71] = soln[10]
        algebraic[72] = soln[11]
        algebraic[73] = soln[12]
        algebraic[74] = soln[13]
        algebraic[75] = soln[14]
        algebraic[76] = soln[15]
        algebraic[77] = soln[16]
    else:
        for (i,t) in enumerate(voi):
            soln = fsolve(residualSN_1, initialGuess1, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6)
            initialGuess1 = soln
            algebraic[61][i] = soln[0]
            algebraic[62][i] = soln[1]
            algebraic[63][i] = soln[2]
            algebraic[64][i] = soln[3]
            algebraic[65][i] = soln[4]
            algebraic[66][i] = soln[5]
            algebraic[67][i] = soln[6]
            algebraic[68][i] = soln[7]
            algebraic[69][i] = soln[8]
            algebraic[70][i] = soln[9]
            algebraic[71][i] = soln[10]
            algebraic[72][i] = soln[11]
            algebraic[73][i] = soln[12]
            algebraic[74][i] = soln[13]
            algebraic[75][i] = soln[14]
            algebraic[76][i] = soln[15]
            algebraic[77][i] = soln[16]

def residualSN_1(algebraicCandidate, algebraic, voi, constants, rates, states):
    resid = array([0.0] * 17)
    algebraic[61] = algebraicCandidate[0]
    algebraic[62] = algebraicCandidate[1]
    algebraic[63] = algebraicCandidate[2]
    algebraic[64] = algebraicCandidate[3]
    algebraic[65] = algebraicCandidate[4]
    algebraic[66] = algebraicCandidate[5]
    algebraic[67] = algebraicCandidate[6]
    algebraic[68] = algebraicCandidate[7]
    algebraic[69] = algebraicCandidate[8]
    algebraic[70] = algebraicCandidate[9]
    algebraic[71] = algebraicCandidate[10]
    algebraic[72] = algebraicCandidate[11]
    algebraic[73] = algebraicCandidate[12]
    algebraic[74] = algebraicCandidate[13]
    algebraic[75] = algebraicCandidate[14]
    algebraic[76] = algebraicCandidate[15]
    algebraic[77] = algebraicCandidate[16]
    resid[0] = (algebraic[61]-(algebraic[64]*algebraic[65])/constants[4])
    resid[1] = (algebraic[62]-(algebraic[61]*algebraic[66])/constants[5])
    resid[2] = (algebraic[63]-(algebraic[65]*algebraic[66])/constants[6])
    resid[3] = (algebraic[64]-((constants[1]-algebraic[61])-algebraic[62]))
    resid[4] = (algebraic[65]-(((states[0]-algebraic[61])-algebraic[62])-algebraic[63]))
    resid[5] = (algebraic[66]-((constants[3]-algebraic[62])-algebraic[63]))
    resid[6] = (algebraic[69]-(constants[32]*constants[36])/(constants[36]+algebraic[67]+algebraic[77]))
    resid[7] = (algebraic[70]-(algebraic[67]/constants[35])*algebraic[67]*(1.00000+algebraic[69]/constants[36]))
    resid[8] = (algebraic[71]-algebraic[67]*(1.00000+algebraic[69]/constants[36]))
    resid[9] = (algebraic[72]-(algebraic[77]/constants[35])*algebraic[77]*(1.00000+algebraic[69]/constants[36]))
    resid[10] = (algebraic[73]-algebraic[77]*(1.00000+algebraic[69]/constants[36]))
    resid[11] = (algebraic[74]-(constants[33]/algebraic[68])*algebraic[70])
    resid[12] = (algebraic[75]-(constants[33]/algebraic[68])*algebraic[72])
    resid[13] = (algebraic[76]-((constants[33]*constants[34])/constants[35]+(constants[33]*algebraic[68])/constants[35]+(power(algebraic[68], 2.00000))/constants[35]))
    resid[14] = (algebraic[68]-((states[6]-(algebraic[74]+2.00000*algebraic[70]+2.00000*algebraic[71]))-(algebraic[75]+2.00000*algebraic[72]+2.00000*algebraic[73])))
    resid[15] = (constants[112]-(2.00000*constants[30]*(power(algebraic[68], 2.00000))-algebraic[67]*(1.00000+algebraic[69]/constants[36])*(algebraic[76]*algebraic[67]+power(algebraic[68], 2.00000))))
    resid[16] = (constants[113]-(2.00000*constants[31]*(power(algebraic[68], 2.00000))-algebraic[77]*(1.00000+algebraic[69]/constants[36])*(algebraic[76]*algebraic[77]+power(algebraic[68], 2.00000))))
    return resid

initialGuess2 = None
def rootfind_2(voi, constants, rates, states, algebraic):
    """Calculate values of algebraic variables for DAE"""
    from scipy.optimize import fsolve
    global initialGuess2
    if initialGuess2 is None: initialGuess2 = ones(3)*0.1
    if not iterable(voi):
        soln = fsolve(residualSN_2, initialGuess2, args=(algebraic, voi, constants, rates, states), xtol=1E-6)
        initialGuess2 = soln
        algebraic[41] = soln[0]
        algebraic[42] = soln[1]
        algebraic[43] = soln[2]
    else:
        for (i,t) in enumerate(voi):
            soln = fsolve(residualSN_2, initialGuess2, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6)
            initialGuess2 = soln
            algebraic[41][i] = soln[0]
            algebraic[42][i] = soln[1]
            algebraic[43][i] = soln[2]

def residualSN_2(algebraicCandidate, algebraic, voi, constants, rates, states):
    resid = array([0.0] * 3)
    algebraic[41] = algebraicCandidate[0]
    algebraic[42] = algebraicCandidate[1]
    algebraic[43] = algebraicCandidate[2]
    resid[0] = (algebraic[43]-(algebraic[41]*algebraic[42])/constants[27])
    resid[1] = (algebraic[41]-(states[3]-algebraic[43]))
    resid[2] = (algebraic[42]-(constants[14]-algebraic[43]))
    return resid

initialGuess3 = None
def rootfind_3(voi, constants, rates, states, algebraic):
    """Calculate values of algebraic variables for DAE"""
    from scipy.optimize import fsolve
    global initialGuess3
    if initialGuess3 is None: initialGuess3 = ones(2)*0.1
    if not iterable(voi):
        soln = fsolve(residualSN_3, initialGuess3, args=(algebraic, voi, constants, rates, states), xtol=1E-6)
        initialGuess3 = soln
        algebraic[49] = soln[0]
        algebraic[50] = soln[1]
    else:
        for (i,t) in enumerate(voi):
            soln = fsolve(residualSN_3, initialGuess3, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6)
            initialGuess3 = soln
            algebraic[49][i] = soln[0]
            algebraic[50][i] = soln[1]

def residualSN_3(algebraicCandidate, algebraic, voi, constants, rates, states):
    resid = array([0.0] * 2)
    algebraic[49] = algebraicCandidate[0]
    algebraic[50] = algebraicCandidate[1]
    resid[0] = (algebraic[50]-(algebraic[49]*algebraic[42])/constants[28])
    resid[1] = (algebraic[49]-(constants[18]-algebraic[50]))
    return resid

initialGuess4 = None
def rootfind_4(voi, constants, rates, states, algebraic):
    """Calculate values of algebraic variables for DAE"""
    from scipy.optimize import fsolve
    global initialGuess4
    if initialGuess4 is None: initialGuess4 = ones(3)*0.1
    if not iterable(voi):
        soln = fsolve(residualSN_4, initialGuess4, args=(algebraic, voi, constants, rates, states), xtol=1E-6)
        initialGuess4 = soln
        algebraic[54] = soln[0]
        algebraic[55] = soln[1]
        algebraic[56] = soln[2]
    else:
        for (i,t) in enumerate(voi):
            soln = fsolve(residualSN_4, initialGuess4, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6)
            initialGuess4 = soln
            algebraic[54][i] = soln[0]
            algebraic[55][i] = soln[1]
            algebraic[56][i] = soln[2]

def residualSN_4(algebraicCandidate, algebraic, voi, constants, rates, states):
    resid = array([0.0] * 3)
    algebraic[54] = algebraicCandidate[0]
    algebraic[55] = algebraicCandidate[1]
    algebraic[56] = algebraicCandidate[2]
    resid[0] = (algebraic[54]-(algebraic[55]*algebraic[56])/constants[48])
    resid[1] = (algebraic[55]-(states[8]-algebraic[54]))
    resid[2] = (algebraic[56]-(constants[38]-algebraic[54]))
    return resid

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)