# Size of variable arrays:
sizeAlgebraic = 64
sizeStates = 36
sizeConstants = 125
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 (s)"
    legend_states[0] = "NAn in component NAn (mM)"
    legend_algebraic[0] = "Vn_leak_Na in component Vn_leak_Na (mM_per_s)"
    legend_algebraic[3] = "Vn_pump in component Vn_pump (mM_per_s)"
    legend_algebraic[39] = "Vn_stim in component Vn_stim (mM_per_s)"
    legend_states[1] = "GLCn in component GLCn (mM)"
    legend_algebraic[4] = "V_en_GLC in component V_en_GLC (mM_per_s)"
    legend_algebraic[5] = "Vn_hk in component Vn_hk (mM_per_s)"
    legend_states[2] = "G6Pn in component G6Pn (mM)"
    legend_algebraic[6] = "Vn_pgi in component Vn_pgi (mM_per_s)"
    legend_states[3] = "F6Pn in component F6Pn (mM)"
    legend_algebraic[7] = "Vn_pfk in component Vn_pfk (mM_per_s)"
    legend_states[4] = "GAPn in component GAPn (mM)"
    legend_algebraic[47] = "Vn_pgk in component Vn_pgk (mM_per_s)"
    legend_states[5] = "PEPn in component PEPn (mM)"
    legend_algebraic[49] = "Vn_pk in component Vn_pk (mM_per_s)"
    legend_states[6] = "PYRn in component PYRn (mM)"
    legend_algebraic[40] = "Vn_ldh in component Vn_ldh (mM_per_s)"
    legend_algebraic[50] = "Vn_mito in component Vn_mito (mM_per_s)"
    legend_states[7] = "LACn in component LACn (mM)"
    legend_algebraic[8] = "Vne_LAC in component Vne_LAC (mM_per_s)"
    legend_states[8] = "NADHn in component NADHn (mM)"
    legend_states[9] = "ATPn in component ATPn (mM)"
    legend_constants[0] = "nOP in component model_parameters (dimensionless)"
    legend_algebraic[9] = "Vn_ATPase in component Vn_ATPase (mM_per_s)"
    legend_algebraic[51] = "Vn_ck in component Vn_ck (mM_per_s)"
    legend_algebraic[55] = "dAMP_dATPn in component dAMP_dATPn (dimensionless)"
    legend_states[10] = "PCrn in component PCrn (mM)"
    legend_states[11] = "O2n in component O2n (mM)"
    legend_constants[1] = "NAero in component model_parameters (dimensionless)"
    legend_algebraic[10] = "Vcn_O2 in component Vcn_O2 (mM_per_s)"
    legend_states[12] = "GLUn in component GLUn (mM)"
    legend_constants[2] = "Rng in component model_parameters (dimensionless)"
    legend_algebraic[28] = "Vg_gs in component Vg_gs (mM_per_s)"
    legend_algebraic[42] = "Vn_stim_GLU in component Vn_stim_GLU (mM_per_s)"
    legend_states[13] = "NAg in component NAg (mM)"
    legend_algebraic[11] = "Vg_leak_Na in component Vg_leak_Na (mM_per_s)"
    legend_algebraic[12] = "Vg_pump in component Vg_pump (mM_per_s)"
    legend_algebraic[29] = "Veg_GLU in component Veg_GLU (mM_per_s)"
    legend_states[14] = "GLCg in component GLCg (mM)"
    legend_algebraic[14] = "Vcg_GLC in component Vcg_GLC (mM_per_s)"
    legend_algebraic[13] = "Veg_GLC in component Veg_GLC (mM_per_s)"
    legend_algebraic[15] = "Vg_hk in component Vg_hk (mM_per_s)"
    legend_states[15] = "G6Pg in component G6Pg (mM)"
    legend_algebraic[16] = "Vg_pgi in component Vg_pgi (mM_per_s)"
    legend_algebraic[18] = "Vg_glys in component Vg_glys (mM_per_s)"
    legend_algebraic[24] = "Vg_glyp in component Vg_glyp (mM_per_s)"
    legend_states[16] = "F6Pg in component F6Pg (mM)"
    legend_algebraic[17] = "Vg_pfk in component Vg_pfk (mM_per_s)"
    legend_states[17] = "GAPg in component GAPg (mM)"
    legend_algebraic[56] = "Vg_pgk in component Vg_pgk (mM_per_s)"
    legend_states[18] = "PEPg in component PEPg (mM)"
    legend_algebraic[58] = "Vg_pk in component Vg_pk (mM_per_s)"
    legend_states[19] = "PYRg in component PYRg (mM)"
    legend_algebraic[43] = "Vg_ldh in component Vg_ldh (mM_per_s)"
    legend_algebraic[59] = "Vg_mito in component Vg_mito (mM_per_s)"
    legend_states[20] = "LACg in component LACg (mM)"
    legend_algebraic[19] = "Vge_LAC in component Vge_LAC (mM_per_s)"
    legend_algebraic[21] = "Vgc_LAC in component Vgc_LAC (mM_per_s)"
    legend_states[21] = "NADHg in component NADHg (mM)"
    legend_states[22] = "ATPg in component ATPg (mM)"
    legend_algebraic[23] = "Vg_ATPase in component Vg_ATPase (mM_per_s)"
    legend_algebraic[60] = "Vg_ck in component Vg_ck (mM_per_s)"
    legend_algebraic[63] = "dAMP_dATPg in component dAMP_dATPg (dimensionless)"
    legend_states[23] = "PCrg in component PCrg (mM)"
    legend_states[24] = "O2g in component O2g (mM)"
    legend_algebraic[25] = "Vcg_O2 in component Vcg_O2 (mM_per_s)"
    legend_states[25] = "GLYg in component GLYg (mM)"
    legend_states[26] = "GLUg in component GLUg (mM)"
    legend_states[27] = "GLCe in component GLCe (mM)"
    legend_constants[3] = "Reg in component model_parameters (dimensionless)"
    legend_constants[4] = "Ren in component model_parameters (dimensionless)"
    legend_algebraic[26] = "Vce_GLC in component Vce_GLC (mM_per_s)"
    legend_states[28] = "LACe in component LACe (mM)"
    legend_algebraic[27] = "Vec_LAC in component Vec_LAC (mM_per_s)"
    legend_states[29] = "GLUe in component GLUe (mM)"
    legend_states[30] = "O2c in component O2c (mM)"
    legend_constants[5] = "Rcn in component model_parameters (dimensionless)"
    legend_constants[6] = "Rcg in component model_parameters (dimensionless)"
    legend_algebraic[33] = "Vc_O2 in component Vc_O2 (mM_per_s)"
    legend_states[31] = "GLCc in component GLCc (mM)"
    legend_constants[7] = "Rce in component model_parameters (dimensionless)"
    legend_algebraic[34] = "Vc_GLC in component Vc_GLC (mM_per_s)"
    legend_states[32] = "LACc in component LACc (mM)"
    legend_algebraic[35] = "Vc_LAC in component Vc_LAC (mM_per_s)"
    legend_states[33] = "CO2c in component CO2c (mM)"
    legend_algebraic[52] = "Vnc_CO2 in component Vnc_CO2 (mM_per_s)"
    legend_algebraic[32] = "Vc_CO2 in component Vc_CO2 (mM_per_s)"
    legend_algebraic[61] = "Vgc_CO2 in component Vgc_CO2 (mM_per_s)"
    legend_states[34] = "Vv in component Vv (dimensionless)"
    legend_algebraic[30] = "Fin_t in component Fin_t (per_s)"
    legend_algebraic[36] = "Fout_t in component Fout_t (per_s)"
    legend_states[35] = "dHb in component dHb (mM)"
    legend_constants[8] = "O2a in component model_parameters (mM)"
    legend_constants[9] = "gn_NA in component Vn_leak_Na (mS_per_cm2)"
    legend_constants[10] = "Sm_n in component model_parameters (per_cm)"
    legend_constants[11] = "Vm in component model_parameters (mV)"
    legend_constants[12] = "Vn in component model_parameters (dimensionless)"
    legend_constants[13] = "RT in component model_parameters (mV_C_per_mol)"
    legend_constants[14] = "F in component model_parameters (C_per_mole)"
    legend_constants[15] = "NAe in component model_parameters (mM)"
    legend_constants[16] = "kpump in component model_parameters (cm_per_mM_per_s)"
    legend_constants[17] = "Km_pump in component model_parameters (mM)"
    legend_algebraic[37] = "v_stim in component v_stim (mM_per_s)"
    legend_constants[18] = "Km_en_GLC in component V_en_GLC (mM)"
    legend_constants[19] = "Vm_en_GLC in component V_en_GLC (mM_per_s)"
    legend_constants[20] = "Vmax_n_hk in component Vn_hk (mM_per_s)"
    legend_constants[21] = "Km_GLC in component model_parameters (mM)"
    legend_constants[22] = "G6P_inh_hk in component model_parameters (mM)"
    legend_constants[23] = "aG6P_inh_hk in component model_parameters (dimensionless)"
    legend_constants[24] = "Vmaxf_n_pgi in component Vn_pgi (mM_per_s)"
    legend_constants[25] = "Vmaxr_n_pgi in component Vn_pgi (mM_per_s)"
    legend_constants[26] = "Km_G6P in component model_parameters (mM)"
    legend_constants[27] = "Km_F6P_pgi in component model_parameters (mM)"
    legend_constants[28] = "kn_pfk in component Vn_pfk (per_s)"
    legend_constants[29] = "Km_F6P_pfk in component model_parameters (mM)"
    legend_constants[30] = "Ki_ATP in component model_parameters (mM)"
    legend_constants[31] = "nH in component model_parameters (dimensionless)"
    legend_constants[32] = "kn_pgk in component Vn_pgk (per_mM_per_s)"
    legend_algebraic[46] = "ADPn in component ADPn (mM)"
    legend_algebraic[38] = "NADn in component NADn (mM)"
    legend_constants[33] = "kn_pk in component Vn_pk (per_mM_per_s)"
    legend_constants[34] = "kfn_ldh in component Vn_ldh (per_mM_per_s)"
    legend_constants[35] = "krn_ldh in component Vn_ldh (per_mM_per_s)"
    legend_constants[36] = "Vmax_n_mito in component Vn_mito (mM_per_s)"
    legend_constants[37] = "Km_O2 in component model_parameters (mM)"
    legend_constants[38] = "Km_ADP in component model_parameters (mM)"
    legend_constants[39] = "Km_PYR in component model_parameters (mM)"
    legend_constants[40] = "rATP_mito in component model_parameters (dimensionless)"
    legend_constants[41] = "aATP_mito in component model_parameters (dimensionless)"
    legend_constants[42] = "Vmax_ne_LAC in component Vne_LAC (mM_per_s)"
    legend_constants[43] = "Km_ne_LAC in component Vne_LAC (mM)"
    legend_constants[44] = "Vmax_n_ATPase in component Vn_ATPase (mM_per_s)"
    legend_constants[45] = "krn_ck in component Vn_ck (per_mM_per_s)"
    legend_constants[46] = "kfn_ck in component Vn_ck (per_mM_per_s)"
    legend_algebraic[44] = "CRn in component CRn (mM)"
    legend_constants[47] = "nh_O2 in component Vcn_O2 (dimensionless)"
    legend_constants[48] = "PScapn in component Vcn_O2 (per_s)"
    legend_constants[49] = "Ko2 in component model_parameters (mM)"
    legend_constants[50] = "HbOP in component model_parameters (mM)"
    legend_constants[51] = "gg_NA in component Vg_leak_Na (mS_per_cm2)"
    legend_constants[52] = "Sm_g in component model_parameters (per_cm)"
    legend_constants[53] = "Vg in component model_parameters (dimensionless)"
    legend_constants[54] = "Km_eg_GLC in component Veg_GLC (mM)"
    legend_constants[55] = "Vm_eg_GLC in component Veg_GLC (mM_per_s)"
    legend_constants[56] = "KO1 in component model_parameters (dimensionless)"
    legend_constants[57] = "Km_cg_GLC in component Vcg_GLC (mM)"
    legend_constants[58] = "Vm_cg_GLC in component Vcg_GLC (mM_per_s)"
    legend_constants[59] = "Vmax_g_hk in component Vg_hk (mM_per_s)"
    legend_constants[60] = "Vmaxf_g_pgi in component Vg_pgi (mM_per_s)"
    legend_constants[61] = "Vmaxr_g_pgi in component Vg_pgi (mM_per_s)"
    legend_constants[62] = "kg_pfk in component Vg_pfk (per_s)"
    legend_constants[63] = "kg_pgk in component Vg_pgk (per_mM_per_s)"
    legend_algebraic[54] = "ADPg in component ADPg (mM)"
    legend_algebraic[41] = "NADg in component NADg (mM)"
    legend_constants[64] = "kg_pk in component Vg_pk (per_mM_per_s)"
    legend_constants[65] = "kfg_ldh in component Vg_ldh (per_mM_per_s)"
    legend_constants[66] = "krg_ldh in component Vg_ldh (per_mM_per_s)"
    legend_constants[67] = "Vmax_g_mito in component Vg_mito (mM_per_s)"
    legend_constants[68] = "Vmax_ge_LAC in component Vge_LAC (mM_per_s)"
    legend_constants[69] = "Km_ge_LAC in component Vge_LAC (mM)"
    legend_constants[70] = "Vmax_gc_LAC in component Vgc_LAC (mM_per_s)"
    legend_constants[71] = "Km_gc_LAC in component Vgc_LAC (mM)"
    legend_constants[72] = "Vmax_g_ATPase in component Vg_ATPase (mM_per_s)"
    legend_constants[73] = "krg_ck in component Vg_ck (per_mM_per_s)"
    legend_constants[74] = "kfg_ck in component Vg_ck (per_mM_per_s)"
    legend_algebraic[45] = "CRg in component CRg (mM)"
    legend_constants[75] = "PScapg in component Vcg_O2 (per_s)"
    legend_constants[76] = "nh_O2 in component model_parameters (dimensionless)"
    legend_constants[77] = "Vc in component model_parameters (dimensionless)"
    legend_constants[78] = "GLCa in component model_parameters (mM)"
    legend_constants[79] = "Km_ce_GLC in component Vce_GLC (mM)"
    legend_constants[80] = "Vm_ce_GLC in component Vce_GLC (mM_per_s)"
    legend_constants[81] = "LACa in component model_parameters (mM)"
    legend_constants[82] = "Km_ec_LAC in component Vec_LAC (mM)"
    legend_constants[83] = "Vm_ec_LAC in component Vec_LAC (mM_per_s)"
    legend_constants[84] = "R_GLU_NA in component model_parameters (dimensionless)"
    legend_constants[85] = "Km_GLU in component model_parameters (mM)"
    legend_constants[86] = "KO2 in component model_parameters (dimensionless)"
    legend_constants[87] = "Vmax_g_gs in component Vg_gs (mM_per_s)"
    legend_constants[88] = "Km_ATP in component model_parameters (mM)"
    legend_constants[89] = "Vmax_eg_GLU in component Veg_GLU (mM_per_s)"
    legend_constants[90] = "CO2a in component model_parameters (mM)"
    legend_constants[91] = "Vmax_glys in component Vg_glys (mM_per_s)"
    legend_constants[92] = "Km_G6P_glys in component Vg_glys (mM)"
    legend_constants[93] = "GLY_inh in component model_parameters (mM)"
    legend_constants[94] = "aGLY_inh in component model_parameters (dimensionless)"
    legend_constants[95] = "Vmax_glyp in component Vg_glyp (mM_per_s)"
    legend_constants[96] = "Km_GLY in component Vg_glyp (mM)"
    legend_algebraic[22] = "deltaVt_GLY in component Vg_glyp (dimensionless)"
    legend_algebraic[20] = "unitstepSB2 in component unitstepSB2 (dimensionless)"
    legend_constants[97] = "stim in component model_parameters (dimensionless)"
    legend_constants[98] = "to in component model_parameters (s)"
    legend_constants[99] = "to_GLY in component model_parameters (s)"
    legend_constants[100] = "tend_GLY in component model_parameters (s)"
    legend_constants[101] = "sr_GLY in component model_parameters (dimensionless)"
    legend_constants[102] = "t1 in component model_parameters (s)"
    legend_constants[103] = "delta_GLY in component model_parameters (dimensionless)"
    legend_constants[104] = "KO3 in component model_parameters (dimensionless)"
    legend_constants[105] = "CBF0 in component Fin_t (per_s)"
    legend_constants[106] = "tend in component model_parameters (s)"
    legend_constants[107] = "sr in component model_parameters (dimensionless)"
    legend_constants[108] = "deltaf in component model_parameters (dimensionless)"
    legend_constants[109] = "CBF0 in component model_parameters (per_s)"
    legend_constants[110] = "Vv0 in component model_parameters (dimensionless)"
    legend_constants[111] = "tv in component model_parameters (s)"
    legend_constants[112] = "NADH_n_tot in component NADn (mM)"
    legend_constants[113] = "NADH_g_tot in component NADg (mM)"
    legend_constants[114] = "PCrn_tot in component CRn (mM)"
    legend_constants[115] = "PCrg_tot in component CRg (mM)"
    legend_constants[116] = "ATPtot in component model_parameters (mM)"
    legend_constants[117] = "qak in component model_parameters (dimensionless)"
    legend_algebraic[53] = "u_n in component u_n (dimensionless)"
    legend_algebraic[62] = "u_g in component u_g (dimensionless)"
    legend_algebraic[48] = "AMPn in component AMPn (mM)"
    legend_algebraic[57] = "AMPg in component AMPg (mM)"
    legend_algebraic[1] = "BOLD in component BOLD (dimensionless)"
    legend_constants[118] = "k1 in component model_parameters (dimensionless)"
    legend_constants[119] = "k2 in component model_parameters (dimensionless)"
    legend_constants[120] = "k3 in component model_parameters (dimensionless)"
    legend_constants[121] = "dHb0 in component model_parameters (mM)"
    legend_algebraic[31] = "unitpulseSB in component v_stim (dimensionless)"
    legend_constants[122] = "t_n_stim in component model_parameters (s)"
    legend_constants[123] = "v1_n in component model_parameters (mM_per_s)"
    legend_constants[124] = "v2_n in component model_parameters (mM_per_s)"
    legend_algebraic[2] = "unitstepSB in component unitstepSB (dimensionless)"
    legend_rates[0] = "d/dt NAn in component NAn (mM)"
    legend_rates[1] = "d/dt GLCn in component GLCn (mM)"
    legend_rates[2] = "d/dt G6Pn in component G6Pn (mM)"
    legend_rates[3] = "d/dt F6Pn in component F6Pn (mM)"
    legend_rates[4] = "d/dt GAPn in component GAPn (mM)"
    legend_rates[5] = "d/dt PEPn in component PEPn (mM)"
    legend_rates[6] = "d/dt PYRn in component PYRn (mM)"
    legend_rates[7] = "d/dt LACn in component LACn (mM)"
    legend_rates[8] = "d/dt NADHn in component NADHn (mM)"
    legend_rates[9] = "d/dt ATPn in component ATPn (mM)"
    legend_rates[10] = "d/dt PCrn in component PCrn (mM)"
    legend_rates[11] = "d/dt O2n in component O2n (mM)"
    legend_rates[12] = "d/dt GLUn in component GLUn (mM)"
    legend_rates[13] = "d/dt NAg in component NAg (mM)"
    legend_rates[14] = "d/dt GLCg in component GLCg (mM)"
    legend_rates[15] = "d/dt G6Pg in component G6Pg (mM)"
    legend_rates[16] = "d/dt F6Pg in component F6Pg (mM)"
    legend_rates[17] = "d/dt GAPg in component GAPg (mM)"
    legend_rates[18] = "d/dt PEPg in component PEPg (mM)"
    legend_rates[19] = "d/dt PYRg in component PYRg (mM)"
    legend_rates[20] = "d/dt LACg in component LACg (mM)"
    legend_rates[21] = "d/dt NADHg in component NADHg (mM)"
    legend_rates[22] = "d/dt ATPg in component ATPg (mM)"
    legend_rates[23] = "d/dt PCrg in component PCrg (mM)"
    legend_rates[24] = "d/dt O2g in component O2g (mM)"
    legend_rates[25] = "d/dt GLYg in component GLYg (mM)"
    legend_rates[26] = "d/dt GLUg in component GLUg (mM)"
    legend_rates[27] = "d/dt GLCe in component GLCe (mM)"
    legend_rates[28] = "d/dt LACe in component LACe (mM)"
    legend_rates[29] = "d/dt GLUe in component GLUe (mM)"
    legend_rates[30] = "d/dt O2c in component O2c (mM)"
    legend_rates[31] = "d/dt GLCc in component GLCc (mM)"
    legend_rates[32] = "d/dt LACc in component LACc (mM)"
    legend_rates[33] = "d/dt CO2c in component CO2c (mM)"
    legend_rates[34] = "d/dt Vv in component Vv (dimensionless)"
    legend_rates[35] = "d/dt dHb in component dHb (mM)"
    return (legend_states, legend_algebraic, legend_voi, legend_constants)

def initConsts():
    constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
    states[0] = 15.533
    states[1] = 0.2633
    states[2] = 0.7275
    states[3] = 0.1091
    states[4] = 0.0418
    states[5] = 0.0037
    states[6] = 0.0388
    states[7] = 0.3856
    states[8] = 0.0319
    states[9] = 2.2592
    constants[0] = 15.0
    states[10] = 4.2529
    states[11] = 0.0975
    constants[1] = 3.0
    states[12] = 3.0
    constants[2] = 1.8
    states[13] = 13.36
    states[14] = 0.1656
    states[15] = 0.7326
    states[16] = 0.1116
    states[17] = 0.0698
    states[18] = 0.0254
    states[19] = 0.1711
    states[20] = 0.4651
    states[21] = 0.0445
    states[22] = 2.24
    states[23] = 4.6817
    states[24] = 0.1589
    states[25] = 2.5
    states[26] = 0.0
    states[27] = 0.3339
    constants[3] = 0.8
    constants[4] = 0.4444444444444444
    states[28] = 0.3986
    states[29] = 0.0
    states[30] = 7.4201
    constants[5] = 0.01222
    constants[6] = 0.022
    states[31] = 4.6401
    constants[7] = 0.0275
    states[32] = 0.3251
    states[33] = 2.12
    states[34] = 0.0237
    states[35] = 0.0218
    constants[8] = 8.34
    constants[9] = 0.0039
    constants[10] = 40500
    constants[11] = -70
    constants[12] = 0.45
    constants[13] = 2577340
    constants[14] = 96500
    constants[15] = 150.0
    constants[16] = 3.17e-7
    constants[17] = 0.4243
    constants[18] = 5.32
    constants[19] = 0.50417
    constants[20] = 0.0513
    constants[21] = 0.105
    constants[22] = 0.6
    constants[23] = 20.0
    constants[24] = 0.5
    constants[25] = 0.45
    constants[26] = 0.5
    constants[27] = 0.06
    constants[28] = 0.55783
    constants[29] = 0.18
    constants[30] = 0.7595
    constants[31] = 4.0
    constants[32] = 0.4287
    constants[33] = 28.6
    constants[34] = 5.30
    constants[35] = 0.1046
    constants[36] = 0.05557
    constants[37] = 0.0029658
    constants[38] = 0.00107
    constants[39] = 0.0632
    constants[40] = 20.0
    constants[41] = 5.0
    constants[42] = 0.1978
    constants[43] = 0.09314
    constants[44] = 0.04889
    constants[45] = 0.015
    constants[46] = 0.0524681
    constants[47] = 2.7
    constants[48] = 0.2202
    constants[49] = 0.089733
    constants[50] = 8.6
    constants[51] = 0.00325
    constants[52] = 10500
    constants[53] = 0.25
    constants[54] = 3.53
    constants[55] = 0.038089
    constants[56] = 1.0
    constants[57] = 9.92
    constants[58] = 0.0098394
    constants[59] = 0.050461
    constants[60] = 0.5
    constants[61] = 0.45
    constants[62] = 0.403
    constants[63] = 0.2514
    constants[64] = 2.73
    constants[65] = 6.2613
    constants[66] = 0.54682
    constants[67] = 0.008454
    constants[68] = 0.086124
    constants[69] = 0.22163
    constants[70] = 0.00021856
    constants[71] = 0.12862
    constants[72] = 0.035657
    constants[73] = 0.02073
    constants[74] = 0.0243
    constants[75] = 0.2457
    constants[76] = 2.7
    constants[77] = 0.0055
    constants[78] = 4.8
    constants[79] = 8.4568
    constants[80] = 0.0489
    constants[81] = 0.313
    constants[82] = 0.764818
    constants[83] = 0.0325
    constants[84] = 0.075
    constants[85] = 0.05
    constants[86] = 1
    constants[87] = 0.3
    constants[88] = 0.01532
    constants[89] = 0.0208
    constants[90] = 1.2
    constants[91] = 0.0001528
    constants[92] = 0.5
    constants[93] = 4.2
    constants[94] = 20.0
    constants[95] = 4.922e-5
    constants[96] = 1.0
    constants[97] = 1
    constants[98] = 200
    constants[99] = 83
    constants[100] = 440
    constants[101] = 4
    constants[102] = 2
    constants[103] = 62
    constants[104] = 1
    constants[105] = 0.012
    constants[106] = 300
    constants[107] = 4.59186
    constants[108] = 0.42
    constants[109] = 0.012
    constants[110] = 0.0236
    constants[111] = 35.0
    constants[112] = 0.22
    constants[113] = 0.22
    constants[114] = 5.0
    constants[115] = 5.0
    constants[116] = 2.379
    constants[117] = 0.92
    constants[118] = 2.22
    constants[119] = 0.46
    constants[120] = 0.43
    constants[121] = 0.064
    constants[122] = 2
    constants[123] = 0.041
    constants[124] = 2.55
    return (states, constants)

def computeRates(voi, states, constants):
    rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
    algebraic[4] = constants[19]*(states[27]/(states[27]+constants[18])-states[1]/(states[1]+constants[18]))
    algebraic[5] = constants[20]*states[9]*(states[1]/(states[1]+constants[21]))*(1.00000-1.00000/(1.00000+exp(-constants[23]*(1.00000*(states[2]-constants[22])))))
    rates[1] = algebraic[4]-algebraic[5]
    algebraic[6] = constants[24]*(states[2]/(states[2]+constants[26]))-constants[25]*(states[3]/(states[3]+constants[27]))
    rates[2] = algebraic[5]-algebraic[6]
    algebraic[7] = constants[28]*states[9]*(states[3]/(states[3]+constants[29]))*(power(1.00000+power(states[9]/constants[30], constants[31]), -1.00000))
    rates[3] = algebraic[6]-algebraic[7]
    algebraic[14] = constants[58]*(states[31]/(states[31]+constants[57])-states[14]/(states[14]+constants[57]))
    algebraic[13] = constants[56]*constants[55]*(states[27]/(states[27]+constants[54])-states[14]/(states[14]+constants[54]))
    algebraic[15] = constants[59]*states[22]*(states[14]/(states[14]+constants[21]))*(1.00000-1.00000/(1.00000+exp(-constants[23]*(1.00000*(states[15]-constants[22])))))
    rates[14] = (algebraic[14]+algebraic[13])-algebraic[15]
    algebraic[16] = constants[60]*(states[15]/(states[15]+constants[26]))-constants[61]*(states[16]/(states[16]+constants[27]))
    algebraic[17] = constants[62]*states[22]*(states[16]/(states[16]+constants[29]))*(power(1.00000+power(states[22]/constants[30], constants[31]), -1.00000))
    rates[16] = algebraic[16]-algebraic[17]
    algebraic[18] = constants[91]*(states[15]/(states[15]+constants[92]))*(1.00000-1.00000/(1.00000+exp(-constants[94]*(1.00000*(states[25]-constants[93])))))
    algebraic[20] = custom_piecewise([greater_equal(voi-(constants[100]+constants[98]+constants[99]) , 0.00000), 1.00000 , True, 0.00000])
    algebraic[22] = 1.00000+constants[97]*(constants[103]*constants[104]*(1.00000/(1.00000+exp(1.00000*-constants[101]*(voi-(constants[98]+constants[99])))))*(1.00000-algebraic[20]))
    algebraic[24] = constants[95]*(states[25]/(states[25]+constants[96]))*algebraic[22]
    rates[15] = (algebraic[15]+algebraic[24])-(algebraic[16]+algebraic[18])
    rates[25] = algebraic[18]-algebraic[24]
    algebraic[26] = constants[80]*(states[31]/(states[31]+constants[79])-states[27]/(states[27]+constants[79]))
    rates[27] = algebraic[26]-(algebraic[13]*(1.00000/constants[3])+algebraic[4]*(1.00000/constants[4]))
    algebraic[8] = constants[42]*(states[7]/(states[7]+constants[43])-states[28]/(states[28]+constants[43]))
    algebraic[19] = constants[68]*(states[20]/(states[20]+constants[69])-states[28]/(states[28]+constants[69]))
    algebraic[27] = constants[83]*(states[28]/(states[28]+constants[82])-states[32]/(states[32]+constants[82]))
    rates[28] = (algebraic[8]*(1.00000/constants[4])+algebraic[19]*(1.00000/constants[3]))-algebraic[27]
    algebraic[11] = (constants[52]/constants[53])*(constants[51]/constants[14])*((constants[13]/constants[14])*log(constants[15]/states[13])-constants[11])
    algebraic[12] = (constants[52]/constants[53])*constants[16]*states[22]*states[13]*(power(1.00000+states[22]/constants[17], -1.00000))
    algebraic[29] = constants[89]*(states[29]/(states[29]+constants[85]))
    rates[13] = (algebraic[11]+3.00000*algebraic[29])-3.00000*algebraic[12]
    algebraic[28] = constants[87]*((states[26]/(states[26]+constants[85]))*(states[22]/(states[22]+constants[88])))
    rates[26] = algebraic[29]-algebraic[28]
    algebraic[10] = (constants[48]/constants[12])*(constants[49]*(power(constants[50]/states[30]-1.00000, -1.00000/constants[47]))-states[11])
    algebraic[25] = (constants[75]/constants[53])*(constants[49]*(power(constants[50]/states[30]-1.00000, -1.00000/constants[76]))-states[24])
    algebraic[30] = constants[105]+(constants[97]*constants[105]*constants[108]*(1.00000/(1.00000+exp((1.00000*-constants[107])*(voi-((constants[98]+constants[102])-3.00000)))))-constants[97]*constants[105]*constants[108]*(1.00000/(1.00000+exp((1.00000*-constants[107])*(voi-(constants[98]+constants[106]+constants[102]+3.00000))))))
    algebraic[33] = 2.00000*(algebraic[30]/constants[77])*(constants[8]-states[30])
    rates[30] = algebraic[33]-(algebraic[10]*(1.00000/constants[5])+algebraic[25]*(1.00000/constants[6]))
    algebraic[34] = 2.00000*(algebraic[30]/constants[77])*(constants[78]-states[31])
    rates[31] = algebraic[34]-(algebraic[26]*(1.00000/constants[7])+algebraic[14]*(1.00000/constants[6]))
    algebraic[21] = constants[70]*(states[20]/(states[20]+constants[71])-states[32]/(states[32]+constants[71]))
    algebraic[35] = 2.00000*(algebraic[30]/constants[77])*(constants[81]-states[32])
    rates[32] = algebraic[35]+(algebraic[27]*(1.00000/constants[7])+algebraic[21]*(1.00000/constants[6]))
    algebraic[36] = constants[109]*((power(states[34]/constants[110], 2.00000)+constants[111]*(power(states[34]/constants[110], -0.500000))*(algebraic[30]/constants[110]))/(1.00000+constants[109]*constants[111]*(power(states[34]/constants[110], -0.500000))*(1.00000/constants[110])))
    rates[34] = algebraic[30]-algebraic[36]
    rates[35] = algebraic[30]*(constants[8]-states[30])-algebraic[36]*(states[35]/states[34])
    algebraic[0] = (constants[10]/constants[12])*(constants[9]/constants[14])*((constants[13]/constants[14])*log(constants[15]/states[0])-constants[11])
    algebraic[3] = (constants[10]/constants[12])*constants[16]*states[9]*states[0]*(power(1.00000+states[9]/constants[17], -1.00000))
    algebraic[31] = custom_piecewise([greater_equal(voi , constants[98]) & less_equal(voi , constants[98]+constants[106]), 1.00000 , True, 0.00000])
    algebraic[37] = constants[97]*(constants[123]+constants[124]*((voi-constants[98])/constants[122])*exp(-((voi-constants[98])*(algebraic[31]/constants[122]))))*algebraic[31]
    algebraic[39] = algebraic[37]
    rates[0] = (algebraic[0]+algebraic[39])-3.00000*algebraic[3]
    algebraic[38] = constants[112]-states[8]
    algebraic[40] = constants[34]*states[6]*states[8]-constants[35]*states[7]*algebraic[38]
    rates[7] = algebraic[40]-algebraic[8]
    algebraic[42] = algebraic[39]*constants[84]*constants[86]*(states[12]/(states[12]+constants[85]))
    rates[12] = algebraic[28]*(1.00000/constants[2])-algebraic[42]
    rates[29] = algebraic[42]*(1.00000/constants[4])-algebraic[29]*(1.00000/constants[3])
    algebraic[41] = constants[113]-states[21]
    algebraic[43] = constants[65]*states[19]*states[21]-constants[66]*states[20]*algebraic[41]
    rates[20] = algebraic[43]-(algebraic[19]+algebraic[21])
    algebraic[46] = (states[9]/2.00000)*(-constants[117]+power(power(constants[117], 2.00000)+4.00000*constants[117]*(constants[116]/states[9]-1.00000), 1.0/2))
    algebraic[47] = constants[32]*states[4]*algebraic[46]*(algebraic[38]/states[8])
    rates[4] = 2.00000*algebraic[7]-algebraic[47]
    algebraic[49] = constants[33]*states[5]*algebraic[46]
    rates[5] = algebraic[47]-algebraic[49]
    algebraic[50] = constants[36]*(states[11]/(states[11]+constants[37]))*(algebraic[46]/(algebraic[46]+constants[38]))*(states[6]/(states[6]+constants[39]))*(1.00000-1.00000/(1.00000+exp(-constants[41]*(1.00000*(states[9]/algebraic[46]-1.00000*constants[40])))))
    rates[6] = algebraic[49]-(algebraic[40]+algebraic[50])
    rates[8] = algebraic[47]-(algebraic[40]+algebraic[50])
    rates[11] = algebraic[10]-constants[1]*algebraic[50]
    algebraic[44] = constants[114]-states[10]
    algebraic[51] = constants[46]*states[10]*algebraic[46]-constants[45]*algebraic[44]*states[9]
    rates[10] = -algebraic[51]
    algebraic[9] = constants[44]*(states[9]/(states[9]+0.00100000))
    algebraic[53] = power(constants[117], 2.00000)+4.00000*constants[117]*(constants[116]/states[9]-1.00000)
    algebraic[55] = (constants[117]/2.00000+constants[117]*(constants[116]/(states[9]*(power(algebraic[53], 1.0/2)))))-(1.00000+0.500000*(power(algebraic[53], 1.0/2)))
    rates[9] = ((algebraic[47]+algebraic[49]+constants[0]*algebraic[50]+algebraic[51])-(algebraic[5]+algebraic[7]+algebraic[9]+algebraic[3]))*(power(1.00000-algebraic[55], -1.00000))
    algebraic[54] = (states[22]/2.00000)*(-constants[117]+power(power(constants[117], 2.00000)+4.00000*constants[117]*(constants[116]/states[22]-1.00000), 1.0/2))
    algebraic[56] = constants[63]*states[17]*algebraic[54]*(algebraic[41]/states[21])
    rates[17] = 2.00000*algebraic[17]-algebraic[56]
    algebraic[58] = constants[64]*states[18]*algebraic[54]
    rates[18] = algebraic[56]-algebraic[58]
    algebraic[59] = constants[67]*(states[24]/(states[24]+constants[37]))*(algebraic[54]/(algebraic[54]+constants[38]))*(states[19]/(states[19]+constants[39]))*(1.00000-1.00000/(1.00000+exp(1.00000*-constants[41]*(states[22]/algebraic[54]-1.00000*constants[40]))))
    rates[19] = algebraic[58]-(algebraic[43]+algebraic[59])
    rates[21] = algebraic[56]-(algebraic[43]+algebraic[59])
    rates[24] = algebraic[25]-constants[1]*algebraic[59]
    algebraic[45] = constants[115]-states[23]
    algebraic[60] = constants[74]*states[23]*algebraic[54]-constants[73]*algebraic[45]*states[22]
    rates[23] = -algebraic[60]
    algebraic[52] = 3.00000*algebraic[50]
    algebraic[32] = 2.00000*(algebraic[30]/constants[77])*(states[33]-constants[90])
    algebraic[61] = 3.00000*algebraic[59]
    rates[33] = (algebraic[52]*(1.00000/constants[5])+algebraic[61]*(1.00000/constants[6]))-algebraic[32]
    algebraic[23] = constants[72]*(states[22]/(states[22]+0.00100000))
    algebraic[62] = power(constants[117], 2.00000)+4.00000*constants[117]*(constants[116]/states[22]-1.00000)
    algebraic[63] = (constants[117]/2.00000+constants[117]*(constants[116]/(states[22]*(power(algebraic[62], 1.0/2)))))-(1.00000+0.500000*(power(algebraic[62], 1.0/2)))
    rates[22] = ((algebraic[56]+algebraic[58]+constants[0]*algebraic[59]+algebraic[60])-(algebraic[15]+algebraic[17]+algebraic[23]+algebraic[12]+algebraic[28]))*(power(1.00000-algebraic[63], -1.00000))
    return(rates)

def computeAlgebraic(constants, states, voi):
    algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
    states = array(states)
    voi = array(voi)
    algebraic[4] = constants[19]*(states[27]/(states[27]+constants[18])-states[1]/(states[1]+constants[18]))
    algebraic[5] = constants[20]*states[9]*(states[1]/(states[1]+constants[21]))*(1.00000-1.00000/(1.00000+exp(-constants[23]*(1.00000*(states[2]-constants[22])))))
    algebraic[6] = constants[24]*(states[2]/(states[2]+constants[26]))-constants[25]*(states[3]/(states[3]+constants[27]))
    algebraic[7] = constants[28]*states[9]*(states[3]/(states[3]+constants[29]))*(power(1.00000+power(states[9]/constants[30], constants[31]), -1.00000))
    algebraic[14] = constants[58]*(states[31]/(states[31]+constants[57])-states[14]/(states[14]+constants[57]))
    algebraic[13] = constants[56]*constants[55]*(states[27]/(states[27]+constants[54])-states[14]/(states[14]+constants[54]))
    algebraic[15] = constants[59]*states[22]*(states[14]/(states[14]+constants[21]))*(1.00000-1.00000/(1.00000+exp(-constants[23]*(1.00000*(states[15]-constants[22])))))
    algebraic[16] = constants[60]*(states[15]/(states[15]+constants[26]))-constants[61]*(states[16]/(states[16]+constants[27]))
    algebraic[17] = constants[62]*states[22]*(states[16]/(states[16]+constants[29]))*(power(1.00000+power(states[22]/constants[30], constants[31]), -1.00000))
    algebraic[18] = constants[91]*(states[15]/(states[15]+constants[92]))*(1.00000-1.00000/(1.00000+exp(-constants[94]*(1.00000*(states[25]-constants[93])))))
    algebraic[20] = custom_piecewise([greater_equal(voi-(constants[100]+constants[98]+constants[99]) , 0.00000), 1.00000 , True, 0.00000])
    algebraic[22] = 1.00000+constants[97]*(constants[103]*constants[104]*(1.00000/(1.00000+exp(1.00000*-constants[101]*(voi-(constants[98]+constants[99])))))*(1.00000-algebraic[20]))
    algebraic[24] = constants[95]*(states[25]/(states[25]+constants[96]))*algebraic[22]
    algebraic[26] = constants[80]*(states[31]/(states[31]+constants[79])-states[27]/(states[27]+constants[79]))
    algebraic[8] = constants[42]*(states[7]/(states[7]+constants[43])-states[28]/(states[28]+constants[43]))
    algebraic[19] = constants[68]*(states[20]/(states[20]+constants[69])-states[28]/(states[28]+constants[69]))
    algebraic[27] = constants[83]*(states[28]/(states[28]+constants[82])-states[32]/(states[32]+constants[82]))
    algebraic[11] = (constants[52]/constants[53])*(constants[51]/constants[14])*((constants[13]/constants[14])*log(constants[15]/states[13])-constants[11])
    algebraic[12] = (constants[52]/constants[53])*constants[16]*states[22]*states[13]*(power(1.00000+states[22]/constants[17], -1.00000))
    algebraic[29] = constants[89]*(states[29]/(states[29]+constants[85]))
    algebraic[28] = constants[87]*((states[26]/(states[26]+constants[85]))*(states[22]/(states[22]+constants[88])))
    algebraic[10] = (constants[48]/constants[12])*(constants[49]*(power(constants[50]/states[30]-1.00000, -1.00000/constants[47]))-states[11])
    algebraic[25] = (constants[75]/constants[53])*(constants[49]*(power(constants[50]/states[30]-1.00000, -1.00000/constants[76]))-states[24])
    algebraic[30] = constants[105]+(constants[97]*constants[105]*constants[108]*(1.00000/(1.00000+exp((1.00000*-constants[107])*(voi-((constants[98]+constants[102])-3.00000)))))-constants[97]*constants[105]*constants[108]*(1.00000/(1.00000+exp((1.00000*-constants[107])*(voi-(constants[98]+constants[106]+constants[102]+3.00000))))))
    algebraic[33] = 2.00000*(algebraic[30]/constants[77])*(constants[8]-states[30])
    algebraic[34] = 2.00000*(algebraic[30]/constants[77])*(constants[78]-states[31])
    algebraic[21] = constants[70]*(states[20]/(states[20]+constants[71])-states[32]/(states[32]+constants[71]))
    algebraic[35] = 2.00000*(algebraic[30]/constants[77])*(constants[81]-states[32])
    algebraic[36] = constants[109]*((power(states[34]/constants[110], 2.00000)+constants[111]*(power(states[34]/constants[110], -0.500000))*(algebraic[30]/constants[110]))/(1.00000+constants[109]*constants[111]*(power(states[34]/constants[110], -0.500000))*(1.00000/constants[110])))
    algebraic[0] = (constants[10]/constants[12])*(constants[9]/constants[14])*((constants[13]/constants[14])*log(constants[15]/states[0])-constants[11])
    algebraic[3] = (constants[10]/constants[12])*constants[16]*states[9]*states[0]*(power(1.00000+states[9]/constants[17], -1.00000))
    algebraic[31] = custom_piecewise([greater_equal(voi , constants[98]) & less_equal(voi , constants[98]+constants[106]), 1.00000 , True, 0.00000])
    algebraic[37] = constants[97]*(constants[123]+constants[124]*((voi-constants[98])/constants[122])*exp(-((voi-constants[98])*(algebraic[31]/constants[122]))))*algebraic[31]
    algebraic[39] = algebraic[37]
    algebraic[38] = constants[112]-states[8]
    algebraic[40] = constants[34]*states[6]*states[8]-constants[35]*states[7]*algebraic[38]
    algebraic[42] = algebraic[39]*constants[84]*constants[86]*(states[12]/(states[12]+constants[85]))
    algebraic[41] = constants[113]-states[21]
    algebraic[43] = constants[65]*states[19]*states[21]-constants[66]*states[20]*algebraic[41]
    algebraic[46] = (states[9]/2.00000)*(-constants[117]+power(power(constants[117], 2.00000)+4.00000*constants[117]*(constants[116]/states[9]-1.00000), 1.0/2))
    algebraic[47] = constants[32]*states[4]*algebraic[46]*(algebraic[38]/states[8])
    algebraic[49] = constants[33]*states[5]*algebraic[46]
    algebraic[50] = constants[36]*(states[11]/(states[11]+constants[37]))*(algebraic[46]/(algebraic[46]+constants[38]))*(states[6]/(states[6]+constants[39]))*(1.00000-1.00000/(1.00000+exp(-constants[41]*(1.00000*(states[9]/algebraic[46]-1.00000*constants[40])))))
    algebraic[44] = constants[114]-states[10]
    algebraic[51] = constants[46]*states[10]*algebraic[46]-constants[45]*algebraic[44]*states[9]
    algebraic[9] = constants[44]*(states[9]/(states[9]+0.00100000))
    algebraic[53] = power(constants[117], 2.00000)+4.00000*constants[117]*(constants[116]/states[9]-1.00000)
    algebraic[55] = (constants[117]/2.00000+constants[117]*(constants[116]/(states[9]*(power(algebraic[53], 1.0/2)))))-(1.00000+0.500000*(power(algebraic[53], 1.0/2)))
    algebraic[54] = (states[22]/2.00000)*(-constants[117]+power(power(constants[117], 2.00000)+4.00000*constants[117]*(constants[116]/states[22]-1.00000), 1.0/2))
    algebraic[56] = constants[63]*states[17]*algebraic[54]*(algebraic[41]/states[21])
    algebraic[58] = constants[64]*states[18]*algebraic[54]
    algebraic[59] = constants[67]*(states[24]/(states[24]+constants[37]))*(algebraic[54]/(algebraic[54]+constants[38]))*(states[19]/(states[19]+constants[39]))*(1.00000-1.00000/(1.00000+exp(1.00000*-constants[41]*(states[22]/algebraic[54]-1.00000*constants[40]))))
    algebraic[45] = constants[115]-states[23]
    algebraic[60] = constants[74]*states[23]*algebraic[54]-constants[73]*algebraic[45]*states[22]
    algebraic[52] = 3.00000*algebraic[50]
    algebraic[32] = 2.00000*(algebraic[30]/constants[77])*(states[33]-constants[90])
    algebraic[61] = 3.00000*algebraic[59]
    algebraic[23] = constants[72]*(states[22]/(states[22]+0.00100000))
    algebraic[62] = power(constants[117], 2.00000)+4.00000*constants[117]*(constants[116]/states[22]-1.00000)
    algebraic[63] = (constants[117]/2.00000+constants[117]*(constants[116]/(states[22]*(power(algebraic[62], 1.0/2)))))-(1.00000+0.500000*(power(algebraic[62], 1.0/2)))
    algebraic[1] = constants[110]*((constants[118]+constants[119])*(1.00000-states[35]/constants[121])-(constants[119]+constants[120])*(1.00000-states[34]/constants[110]))
    algebraic[2] = custom_piecewise([greater_equal(voi-(constants[106]+constants[98]) , 0.00000), 1.00000 , True, 0.00000])
    algebraic[48] = constants[116]-(states[9]+algebraic[46])
    algebraic[57] = constants[116]-(states[22]+algebraic[54])
    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)