# Size of variable arrays:
sizeAlgebraic = 37
sizeStates = 20
sizeConstants = 144
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 (hour)"
    legend_algebraic[4] = "V_K in component membrane (nanoM)"
    legend_algebraic[36] = "fr in component membrane (hertz)"
    legend_algebraic[0] = "v_sPC in component membrane (flux)"
    legend_algebraic[24] = "Res in component membrane (gigaohm)"
    legend_algebraic[18] = "V_rest in component membrane (milliV)"
    legend_algebraic[23] = "V_reset in component membrane (milliV)"
    legend_constants[0] = "V_theta in component membrane (milliV)"
    legend_constants[1] = "Cm in component membrane (nanoF)"
    legend_constants[2] = "R in component membrane (joule_per_mole_kelvin)"
    legend_constants[3] = "T in component membrane (kelvin)"
    legend_constants[4] = "T_abs in component membrane (kelvin)"
    legend_constants[5] = "T_room in component membrane (kelvin)"
    legend_constants[6] = "F in component membrane (coulomb_per_mole)"
    legend_algebraic[5] = "alpha in component membrane (dimensionless)"
    legend_algebraic[1] = "beta in component membrane (dimensionless)"
    legend_algebraic[15] = "beta_a in component membrane (dimensionless)"
    legend_algebraic[16] = "c in component membrane (dimensionless)"
    legend_constants[7] = "P_Ca in component membrane (per_milliM)"
    legend_algebraic[2] = "P_K in component membrane (per_milliM)"
    legend_constants[8] = "P_Na in component membrane (per_milliM)"
    legend_constants[9] = "P_Cl in component membrane (per_milliM)"
    legend_constants[10] = "v_PK in component membrane (per_milliM)"
    legend_constants[11] = "npk in component membrane (dimensionless)"
    legend_constants[12] = "K_PK in component membrane (per_nanoM2)"
    legend_constants[13] = "V_R in component membrane (gigaohm)"
    legend_constants[14] = "K_R in component membrane (milliV)"
    legend_constants[15] = "Ca_ex in component membrane (microM)"
    legend_constants[16] = "Na_ex in component membrane (milliM)"
    legend_constants[17] = "K_ex in component membrane (milliM)"
    legend_constants[18] = "Cl_ex in component membrane (milliM)"
    legend_constants[142] = "Na_in in component membrane (milliM)"
    legend_constants[139] = "K_in in component membrane (milliM)"
    legend_constants[141] = "theta_Na in component membrane (dimensionless)"
    legend_constants[138] = "theta_K in component membrane (dimensionless)"
    legend_constants[19] = "E_Na_0 in component membrane (milliV)"
    legend_constants[140] = "E_Na in component membrane (milliV)"
    legend_constants[20] = "g_Na in component membrane (nanoS)"
    legend_algebraic[7] = "E_Ca in component membrane (milliV)"
    legend_algebraic[11] = "g_Ca in component membrane (nanoS)"
    legend_algebraic[9] = "g_K in component membrane (nanoS)"
    legend_constants[21] = "g_K0 in component membrane (nanoS)"
    legend_constants[22] = "v_gk in component membrane (nanoS)"
    legend_constants[23] = "K_gk in component membrane (nanoM)"
    legend_algebraic[12] = "g_KCa in component membrane (nanoS)"
    legend_constants[24] = "v_KCa in component membrane (nanoS)"
    legend_constants[25] = "K_KCa in component membrane (per_nanoM)"
    legend_constants[26] = "nkca in component membrane (dimensionless)"
    legend_algebraic[29] = "g_L in component membrane (nanoS)"
    legend_constants[27] = "E_L_0 in component membrane (milliV)"
    legend_constants[143] = "E_L in component membrane (milliV)"
    legend_constants[28] = "E_K_0 in component membrane (milliV)"
    legend_constants[136] = "E_K in component membrane (milliV)"
    legend_algebraic[30] = "g_ex in component membrane (nanoS)"
    legend_constants[29] = "E_ex in component membrane (milliV)"
    legend_constants[30] = "v_ex1 in component membrane (nanoS)"
    legend_constants[31] = "v_ex2 in component membrane (nanoS)"
    legend_constants[32] = "K_ex1 in component membrane (picoA_2_5)"
    legend_constants[33] = "K_ex2 in component membrane (per_microM)"
    legend_constants[34] = "nex1 in component membrane (dimensionless)"
    legend_constants[35] = "nex2 in component membrane (dimensionless)"
    legend_constants[36] = "g_inhib in component membrane (nanoS)"
    legend_algebraic[25] = "E_inhib in component membrane (milliV)"
    legend_algebraic[13] = "GABA in component membrane (nanoM)"
    legend_constants[37] = "GABA0 in component membrane (nanoM)"
    legend_constants[38] = "v_GABA in component membrane (nanoM)"
    legend_constants[39] = "K_GABA in component membrane (nanoM)"
    legend_algebraic[14] = "Cl_in in component membrane (milliM)"
    legend_constants[40] = "Cl0 in component membrane (milliM)"
    legend_constants[41] = "v_Cl1 in component membrane (milliM)"
    legend_constants[42] = "v_Cl2 in component membrane (milliM)"
    legend_constants[43] = "K_Cl1 in component membrane (nanoM)"
    legend_constants[44] = "K_Cl2 in component membrane (per_nanoM_0_2)"
    legend_constants[45] = "nCl in component membrane (dimensionless)"
    legend_algebraic[26] = "i_Na in component membrane (picoA)"
    legend_algebraic[28] = "i_Na_abs in component membrane (picoA)"
    legend_algebraic[20] = "i_Ca in component membrane (picoA)"
    legend_algebraic[21] = "i_K in component membrane (picoA)"
    legend_algebraic[22] = "i_KCa in component membrane (picoA)"
    legend_algebraic[31] = "i_L in component membrane (picoA)"
    legend_algebraic[32] = "i_ex in component membrane (picoA)"
    legend_algebraic[27] = "i_inhib in component membrane (picoA)"
    legend_algebraic[33] = "I_star in component membrane (picoA)"
    legend_algebraic[34] = "R_star in component membrane (milliohm)"
    legend_algebraic[35] = "tau_m in component membrane (second)"
    legend_algebraic[19] = "theta in component membrane (milliV)"
    legend_algebraic[17] = "psi in component membrane (dimensionless)"
    legend_constants[46] = "k_q in component membrane (per_kelvin)"
    legend_constants[47] = "k_MK in component membrane (nanoM)"
    legend_constants[48] = "V_MK in component membrane (nanoM)"
    legend_constants[49] = "V_b in component membrane (nanoM)"
    legend_constants[50] = "k_b in component membrane (dimensionless)"
    legend_constants[51] = "K_D in component membrane (nanoM)"
    legend_constants[52] = "v_sP0 in component membrane (flux)"
    legend_constants[53] = "C_T in component membrane (flux)"
    legend_constants[54] = "nca in component membrane (dimensionless)"
    legend_constants[55] = "v_Ca in component membrane (nanoS)"
    legend_constants[56] = "K_Ca in component membrane (nanoM_2_2)"
    legend_constants[57] = "K_C in component membrane (nanoM)"
    legend_states[0] = "Ca_in in component Ca_in (microM)"
    legend_states[1] = "BC in component BC (nanoM)"
    legend_states[2] = "MP in component MP (nanoM)"
    legend_states[3] = "VIP in component VIP (nanoM)"
    legend_states[4] = "CC in component CC (nanoM)"
    legend_states[5] = "CB in component CB (nanoM)"
    legend_constants[58] = "beta_IP3 in component Ca_in (dimensionless)"
    legend_states[6] = "Ca_store in component Ca_in (nanoM)"
    legend_algebraic[6] = "v0 in component Ca_in (microM_per_hour)"
    legend_constants[137] = "v1 in component Ca_in (microM_per_hour)"
    legend_algebraic[8] = "v2 in component Ca_in (microM_per_hour)"
    legend_algebraic[10] = "v3 in component Ca_in (microM_per_hour)"
    legend_constants[59] = "v_kk in component Ca_in (per_microM_per_hour)"
    legend_constants[60] = "v_v0 in component Ca_in (microM_per_hour)"
    legend_constants[61] = "kf in component Ca_in (per_hour)"
    legend_algebraic[3] = "k in component Ca_in (per_microM_per_hour)"
    legend_constants[62] = "nkk in component Ca_in (dimensionless)"
    legend_constants[63] = "nv0 in component Ca_in (dimensionless)"
    legend_constants[64] = "m in component Ca_in (dimensionless)"
    legend_constants[65] = "n in component Ca_in (dimensionless)"
    legend_constants[66] = "p in component Ca_in (dimensionless)"
    legend_constants[67] = "v in component Ca_in (dimensionless)"
    legend_constants[68] = "K_kk in component Ca_in (nanoM_0_1)"
    legend_constants[69] = "K_v0 in component Ca_in (nanoM_4_5)"
    legend_constants[70] = "K_2 in component Ca_in (microM)"
    legend_constants[71] = "K_R in component Ca_in (microM)"
    legend_constants[72] = "K_A in component Ca_in (microM)"
    legend_constants[73] = "v_M1 in component Ca_in (microM_per_hour)"
    legend_constants[74] = "v_M2 in component Ca_in (microM_per_hour)"
    legend_constants[75] = "v_M3 in component Ca_in (microM_per_hour)"
    legend_constants[76] = "v_VIP in component VIP (flux)"
    legend_constants[77] = "K_VIP in component VIP (hertz_1_9)"
    legend_constants[78] = "kd_VIP in component VIP (nanoM_0_8_per_hour)"
    legend_constants[79] = "nd_VIP in component VIP (dimensionless)"
    legend_constants[80] = "n_VIP in component VIP (dimensionless)"
    legend_constants[81] = "vmP in component MP (flux)"
    legend_constants[82] = "kdmp in component MP (first_order_rate_constant)"
    legend_constants[83] = "KAP in component MP (nanoM)"
    legend_constants[84] = "KmP in component MP (nanoM)"
    legend_constants[85] = "n in component model_parameters (dimensionless)"
    legend_states[7] = "BN in component BN (nanoM)"
    legend_states[8] = "MC in component MC (nanoM)"
    legend_constants[86] = "vsC in component MC (flux)"
    legend_constants[87] = "vmC in component MC (flux)"
    legend_constants[88] = "kdmc in component MC (first_order_rate_constant)"
    legend_constants[89] = "KAC in component MC (nanoM)"
    legend_constants[90] = "KmC in component MC (nanoM)"
    legend_states[9] = "MB in component MB (nanoM)"
    legend_constants[91] = "vsB in component MB (flux)"
    legend_constants[92] = "vmB in component MB (flux)"
    legend_constants[93] = "kdmb in component MB (first_order_rate_constant)"
    legend_constants[94] = "KIB in component MB (nanoM)"
    legend_constants[95] = "KmB in component MB (nanoM)"
    legend_constants[96] = "m in component model_parameters (dimensionless)"
    legend_states[10] = "PC in component PC (nanoM)"
    legend_constants[97] = "ksP in component model_parameters (first_order_rate_constant)"
    legend_constants[98] = "Kp in component model_parameters (nanoM)"
    legend_constants[99] = "Kdp in component model_parameters (nanoM)"
    legend_constants[100] = "k3 in component model_parameters (second_order_rate_constant)"
    legend_constants[101] = "k4 in component model_parameters (first_order_rate_constant)"
    legend_constants[102] = "kdn in component model_parameters (first_order_rate_constant)"
    legend_constants[103] = "V1P in component model_parameters (flux)"
    legend_constants[104] = "V2P in component model_parameters (flux)"
    legend_states[11] = "PCP in component PCP (nanoM)"
    legend_states[12] = "PCC in component PCC (nanoM)"
    legend_constants[105] = "ksC in component model_parameters (first_order_rate_constant)"
    legend_constants[106] = "kdnc in component model_parameters (first_order_rate_constant)"
    legend_constants[107] = "V1C in component model_parameters (flux)"
    legend_constants[108] = "V2C in component model_parameters (flux)"
    legend_states[13] = "CCP in component CCP (nanoM)"
    legend_constants[109] = "vdPC in component model_parameters (flux)"
    legend_constants[110] = "Kd in component model_parameters (nanoM)"
    legend_constants[111] = "vdCC in component model_parameters (flux)"
    legend_constants[112] = "k1 in component model_parameters (first_order_rate_constant)"
    legend_constants[113] = "k2 in component model_parameters (first_order_rate_constant)"
    legend_constants[114] = "V1PC in component model_parameters (flux)"
    legend_constants[115] = "V2PC in component model_parameters (flux)"
    legend_states[14] = "PCCP in component PCCP (nanoM)"
    legend_states[15] = "PCN in component PCN (nanoM)"
    legend_constants[116] = "k7 in component model_parameters (second_order_rate_constant)"
    legend_constants[117] = "k8 in component model_parameters (first_order_rate_constant)"
    legend_constants[118] = "V3PC in component model_parameters (flux)"
    legend_constants[119] = "V4PC in component model_parameters (flux)"
    legend_states[16] = "PCNP in component PCNP (nanoM)"
    legend_states[17] = "IN in component IN (nanoM)"
    legend_constants[120] = "vdPCC in component model_parameters (flux)"
    legend_constants[121] = "vdPCN in component model_parameters (flux)"
    legend_constants[122] = "ksB in component model_parameters (first_order_rate_constant)"
    legend_constants[123] = "k5 in component model_parameters (first_order_rate_constant)"
    legend_constants[124] = "k6 in component model_parameters (first_order_rate_constant)"
    legend_constants[125] = "V1B in component model_parameters (flux)"
    legend_constants[126] = "V2B in component model_parameters (flux)"
    legend_states[18] = "BCP in component BCP (nanoM)"
    legend_constants[127] = "vdBC in component model_parameters (flux)"
    legend_constants[128] = "V3B in component model_parameters (flux)"
    legend_constants[129] = "V4B in component model_parameters (flux)"
    legend_states[19] = "BNP in component BNP (nanoM)"
    legend_constants[130] = "vdBN in component model_parameters (flux)"
    legend_constants[131] = "vdIN in component model_parameters (flux)"
    legend_constants[132] = "K1CB in component CB (nanoM)"
    legend_constants[133] = "K2CB in component CB (nanoM)"
    legend_constants[134] = "WT in component CB (hour)"
    legend_constants[135] = "vP in component CB (nanoM)"
    legend_rates[0] = "d/dt Ca_in in component Ca_in (microM)"
    legend_rates[6] = "d/dt Ca_store in component Ca_in (nanoM)"
    legend_rates[3] = "d/dt VIP in component VIP (nanoM)"
    legend_rates[2] = "d/dt MP in component MP (nanoM)"
    legend_rates[8] = "d/dt MC in component MC (nanoM)"
    legend_rates[9] = "d/dt MB in component MB (nanoM)"
    legend_rates[10] = "d/dt PC in component PC (nanoM)"
    legend_rates[4] = "d/dt CC in component CC (nanoM)"
    legend_rates[11] = "d/dt PCP in component PCP (nanoM)"
    legend_rates[13] = "d/dt CCP in component CCP (nanoM)"
    legend_rates[12] = "d/dt PCC in component PCC (nanoM)"
    legend_rates[15] = "d/dt PCN in component PCN (nanoM)"
    legend_rates[14] = "d/dt PCCP in component PCCP (nanoM)"
    legend_rates[16] = "d/dt PCNP in component PCNP (nanoM)"
    legend_rates[1] = "d/dt BC in component BC (nanoM)"
    legend_rates[18] = "d/dt BCP in component BCP (nanoM)"
    legend_rates[7] = "d/dt BN in component BN (nanoM)"
    legend_rates[19] = "d/dt BNP in component BNP (nanoM)"
    legend_rates[17] = "d/dt IN in component IN (nanoM)"
    legend_rates[5] = "d/dt CB in component CB (nanoM)"
    return (legend_states, legend_algebraic, legend_voi, legend_constants)

def initConsts():
    constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
    constants[0] = 20.0
    constants[1] = 5.0
    constants[2] = 8.314
    constants[3] = 37
    constants[4] = 273.15
    constants[5] = 22
    constants[6] = 96485
    constants[7] = 0.05
    constants[8] = 0.036
    constants[9] = 0.3
    constants[10] = 1.9
    constants[11] = -2.0
    constants[12] = 1.0
    constants[13] = 0.41
    constants[14] = 34.0
    constants[15] = 5.0
    constants[16] = 145.0
    constants[17] = 1.0
    constants[18] = 114.5
    constants[19] = 45.0
    constants[20] = 36.0
    constants[21] = 9.7
    constants[22] = 10.0
    constants[23] = 10.0
    constants[24] = 3.0
    constants[25] = 0.16
    constants[26] = -1.0
    constants[27] = -29.0
    constants[28] = -97.0
    constants[29] = 0.0
    constants[30] = 105.0
    constants[31] = 4.4
    constants[32] = 574050000.0
    constants[33] = 1.0
    constants[34] = 2.5
    constants[35] = -1.0
    constants[36] = 12.3
    constants[37] = 0.2
    constants[38] = 19.0
    constants[39] = 3.0
    constants[40] = 1.0
    constants[41] = 15.5
    constants[42] = 19.0
    constants[43] = 4.0
    constants[44] = 1.0
    constants[45] = -0.2
    constants[46] = 8.75e-5
    constants[47] = 2.9
    constants[48] = 5.0
    constants[49] = 2.0
    constants[50] = 2.0
    constants[51] = 0.08
    constants[52] = 1.0
    constants[53] = 1.6
    constants[54] = 2.2
    constants[55] = 12.3
    constants[56] = 22.0
    constants[57] = 0.15
    states[0] = 0.1
    states[1] = 2.41
    states[2] = 2.8
    states[3] = 0.0
    states[4] = 12.0
    states[5] = 0.12
    constants[58] = 0.5
    states[6] = 0.1
    constants[59] = 3.3
    constants[60] = 0.09
    constants[61] = 0.001
    constants[62] = 0.1
    constants[63] = 4.5
    constants[64] = 6.0
    constants[65] = 2.2
    constants[66] = 4.2
    constants[67] = 2.0
    constants[68] = 0.02
    constants[69] = 4.5
    constants[70] = 5.0
    constants[71] = 3.0
    constants[72] = 0.67
    constants[73] = 0.0003
    constants[74] = 149.5
    constants[75] = 400.0
    constants[76] = 0.5
    constants[77] = 15.0
    constants[78] = 0.5
    constants[79] = 0.2
    constants[80] = 1.9
    constants[81] = 1.1
    constants[82] = 0.01
    constants[83] = 0.6
    constants[84] = 0.31
    constants[85] = 4.0
    states[7] = 1.94
    states[8] = 2.0
    constants[86] = 1.1
    constants[87] = 1.0
    constants[88] = 0.01
    constants[89] = 0.6
    constants[90] = 0.4
    states[9] = 7.94
    constants[91] = 1.0
    constants[92] = 0.8
    constants[93] = 0.01
    constants[94] = 2.2
    constants[95] = 0.4
    constants[96] = 2.0
    states[10] = 0.4
    constants[97] = 0.6
    constants[98] = 0.1
    constants[99] = 0.1
    constants[100] = 0.4
    constants[101] = 0.2
    constants[102] = 0.01
    constants[103] = 0.4
    constants[104] = 0.3
    states[11] = 0.13
    states[12] = 1.26
    constants[105] = 1.6
    constants[106] = 0.12
    constants[107] = 0.6
    constants[108] = 0.1
    states[13] = 9.0
    constants[109] = 0.7
    constants[110] = 0.3
    constants[111] = 0.7
    constants[112] = 0.45
    constants[113] = 0.2
    constants[114] = 0.4
    constants[115] = 0.1
    states[14] = 0.2
    states[15] = 0.16
    constants[116] = 0.5
    constants[117] = 0.1
    constants[118] = 0.4
    constants[119] = 0.1
    states[16] = 0.091
    states[17] = 0.05
    constants[120] = 0.7
    constants[121] = 0.7
    constants[122] = 0.12
    constants[123] = 0.4
    constants[124] = 0.2
    constants[125] = 0.5
    constants[126] = 0.1
    states[18] = 0.48
    constants[127] = 0.5
    constants[128] = 0.5
    constants[129] = 0.2
    states[19] = 0.32
    constants[130] = 0.6
    constants[131] = 0.8
    constants[132] = 0.01
    constants[133] = 0.01
    constants[134] = 1.0
    constants[135] = 1.0
    constants[136] = (constants[28]*(constants[3]+constants[4]))/(constants[5]+constants[4])
    constants[137] = constants[73]*constants[58]
    constants[138] = exp(constants[136]/(constants[46]*(constants[3]+constants[4])*1000.00))
    constants[139] = constants[17]/constants[138]
    constants[140] = (constants[19]*(constants[3]+constants[4]))/(constants[5]+constants[4])
    constants[141] = exp(constants[140]/(constants[46]*(constants[3]+constants[4])*1000.00))
    constants[142] = constants[16]/constants[141]
    constants[143] = (constants[27]*(constants[3]+constants[4]))/(constants[5]+constants[4])
    return (states, constants)

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

def computeAlgebraic(constants, states, voi):
    algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
    states = array(states)
    voi = array(voi)
    algebraic[0] = constants[52]+(constants[53]*states[5])/(constants[57]+states[5])
    algebraic[1] = states[3]/(states[3]+constants[51])
    algebraic[4] = (constants[48]*states[0])/(constants[47]+states[0])+(constants[49]*algebraic[1])/(constants[50]+algebraic[1])
    algebraic[6] = constants[60]*((power(states[1], constants[63]))/(constants[69]+power(states[1], constants[63])))
    algebraic[8] = constants[74]*((power(states[0], constants[65]))/(power(constants[70], constants[65])+power(states[0], constants[65])))
    algebraic[10] = constants[75]*((power(states[6], constants[64]))/(power(constants[71], constants[64])+power(states[6], constants[64])))*((power(states[0], constants[66]))/(power(constants[72], constants[66])+power(states[0], constants[66])))
    algebraic[3] = constants[59]*((power(states[4], constants[62]))/(constants[68]+power(states[4], constants[62])))
    algebraic[2] = constants[10]*((power(states[1], constants[11]))/(constants[12]+power(states[1], constants[11])))
    algebraic[5] = 4.00000*constants[7]*states[0]*0.00100000+algebraic[2]*constants[139]+constants[8]*constants[142]+constants[9]*constants[18]
    algebraic[13] = constants[37]+(constants[38]*states[3])/(constants[39]+states[3])
    algebraic[14] = constants[40]+constants[41]*(states[2]/(constants[43]+states[2]))+constants[42]*((power(algebraic[13], constants[45]))/(constants[44]+power(algebraic[13], constants[45])))
    algebraic[15] = (algebraic[2]*constants[139]+constants[8]*constants[142]+constants[9]*constants[18])-(algebraic[2]*constants[17]+constants[8]*constants[16]+constants[9]*algebraic[14])
    algebraic[16] = -(algebraic[2]*constants[17]+4.00000*constants[7]*constants[15]*0.00100000+constants[8]*constants[16]+constants[9]*algebraic[14])
    algebraic[17] = (power(power(algebraic[15], 2.00000)-4.00000*algebraic[5]*algebraic[16], 1.0/2)-algebraic[15])/(2.00000*algebraic[5])
    algebraic[18] = ((constants[2]*(constants[3]+constants[4]))/constants[6])*log(algebraic[17])*1000.00
    algebraic[23] = algebraic[18]+4.00000
    algebraic[7] = ((constants[46]*(constants[3]+constants[4]))/2.00000)*log(constants[15]/states[0])*1000.00
    algebraic[11] = constants[55]*((power(states[2], constants[54]))/(constants[56]+power(states[2], constants[54])))
    algebraic[9] = constants[21]+constants[22]*(states[2]/(constants[23]+states[2]))
    algebraic[12] = constants[24]*((power(states[4], constants[26]))/(constants[25]+power(states[4], constants[26])))
    algebraic[24] = constants[13]*(algebraic[18]/(constants[14]+algebraic[18]))
    algebraic[29] = 1.00000/algebraic[24]
    algebraic[26] = constants[20]*(algebraic[18]-constants[140])
    algebraic[28] = power(power(algebraic[26], 2.00000), 1.0/2)
    algebraic[30] = constants[30]*((power(algebraic[28], constants[34]))/(constants[32]+power(algebraic[28], constants[34])))+constants[31]*((power(states[0], constants[35]))/(constants[33]+power(states[0], constants[35])))
    algebraic[25] = -constants[46]*(constants[3]+constants[4])*log(constants[18]/algebraic[14])*1000.00
    algebraic[33] = (constants[20]*constants[140]+algebraic[11]*algebraic[7]+algebraic[9]*constants[136]+algebraic[29]*constants[143]+algebraic[12]*constants[136])-(constants[36]*algebraic[25]+algebraic[30]*constants[29])
    algebraic[34] = 1.00000/((constants[20]+algebraic[9]+algebraic[29]+algebraic[11]+algebraic[12])-(constants[36]+algebraic[30]))
    algebraic[35] = constants[1]*algebraic[34]
    algebraic[19] = algebraic[18]+constants[0]
    algebraic[36] = -1.00000/(algebraic[35]*log((algebraic[19]-algebraic[34]*algebraic[33])/(algebraic[23]-algebraic[34]*algebraic[33])))
    algebraic[20] = algebraic[11]*(algebraic[18]-algebraic[7])
    algebraic[21] = algebraic[9]*(algebraic[18]-constants[136])
    algebraic[22] = algebraic[12]*(algebraic[18]-constants[136])
    algebraic[27] = constants[36]*(algebraic[18]-algebraic[25])
    algebraic[31] = algebraic[29]*(algebraic[18]-constants[143])
    algebraic[32] = algebraic[30]*(algebraic[18]-constants[29])
    return algebraic

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)