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 = 5
sizeStates = 48
sizeConstants = 101
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 (day)"
    legend_algebraic[0] = "FAt in component FAt (litre_per_day)"
    legend_constants[0] = "FA1 in component FAt (litre_per_day)"
    legend_constants[1] = "tA in component FAt (day)"
    legend_constants[2] = "FA2 in component FAt (litre_per_day)"
    legend_algebraic[3] = "FBt in component FBt (litre_per_day)"
    legend_constants[3] = "FB1 in component FBt (litre_per_day)"
    legend_constants[4] = "tB in component FBt (day)"
    legend_constants[5] = "FB2 in component FBt (litre_per_day)"
    legend_algebraic[4] = "F0t in component F0t (litre_per_day)"
    legend_states[0] = "c15 in component c15 (micromole_per_litre)"
    legend_constants[6] = "p0 in component model_parameters (litre_per_day)"
    legend_constants[7] = "V2 in component model_parameters (litre)"
    legend_states[1] = "c015 in component c015 (micromole_per_litre)"
    legend_constants[8] = "p16 in component model_parameters (per_day)"
    legend_states[2] = "c16 in component c16 (micromole_per_litre)"
    legend_constants[9] = "p18 in component model_parameters (per_day)"
    legend_states[3] = "c17 in component c17 (micromole_per_litre)"
    legend_constants[10] = "p15 in component model_parameters (per_day_per_micromole_litre)"
    legend_states[4] = "c19 in component c19 (micromole_per_litre)"
    legend_constants[11] = "p17 in component model_parameters (per_day)"
    legend_constants[12] = "p19 in component model_parameters (per_day_per_micromole_litre)"
    legend_constants[13] = "p20 in component model_parameters (micromole_per_litre)"
    legend_constants[14] = "p29 in component model_parameters (per_day)"
    legend_constants[15] = "pAA15 in component model_parameters (dimensionless)"
    legend_algebraic[1] = "pp in component pp (dimensionless)"
    legend_states[5] = "c016 in component c016 (micromole_per_litre)"
    legend_constants[16] = "pAA16 in component model_parameters (dimensionless)"
    legend_states[6] = "c017 in component c017 (micromole_per_litre)"
    legend_constants[17] = "p22 in component model_parameters (per_day)"
    legend_states[7] = "c18 in component c18 (micromole_per_litre)"
    legend_constants[18] = "p21 in component model_parameters (per_day)"
    legend_constants[19] = "p23 in component model_parameters (per_day)"
    legend_constants[20] = "p1 in component model_parameters (per_day)"
    legend_constants[21] = "s1 in component c17 (dimensionless)"
    legend_states[8] = "c1 in component c1 (micromole_per_litre)"
    legend_constants[22] = "p2 in component model_parameters (per_day)"
    legend_constants[23] = "s2 in component c17 (dimensionless)"
    legend_states[9] = "c2 in component c2 (micromole_per_litre)"
    legend_constants[24] = "p3 in component model_parameters (per_day)"
    legend_constants[25] = "s3 in component c17 (dimensionless)"
    legend_states[10] = "c3 in component c3 (micromole_per_litre)"
    legend_constants[26] = "p4 in component model_parameters (per_day)"
    legend_constants[27] = "s4 in component c17 (dimensionless)"
    legend_states[11] = "c4 in component c4 (micromole_per_litre)"
    legend_constants[28] = "p5 in component model_parameters (per_day)"
    legend_constants[29] = "s5 in component c17 (dimensionless)"
    legend_states[12] = "c5 in component c5 (micromole_per_litre)"
    legend_constants[30] = "p6 in component model_parameters (per_day)"
    legend_constants[31] = "s6 in component c17 (dimensionless)"
    legend_states[13] = "c6 in component c6 (micromole_per_litre)"
    legend_constants[32] = "p7 in component model_parameters (per_day)"
    legend_constants[33] = "s7 in component c17 (dimensionless)"
    legend_states[14] = "c7 in component c7 (micromole_per_litre)"
    legend_constants[34] = "p8 in component model_parameters (per_day)"
    legend_constants[35] = "s8 in component c17 (dimensionless)"
    legend_states[15] = "c8 in component c8 (micromole_per_litre)"
    legend_constants[36] = "p9 in component model_parameters (per_day)"
    legend_constants[37] = "s9 in component c17 (dimensionless)"
    legend_states[16] = "c9 in component c9 (micromole_per_litre)"
    legend_constants[38] = "p14 in component model_parameters (per_day)"
    legend_constants[39] = "s14 in component c17 (dimensionless)"
    legend_states[17] = "c14 in component c14 (micromole_per_litre)"
    legend_constants[40] = "p25 in component model_parameters (per_day_per_micromole_litre)"
    legend_constants[41] = "p26 in component model_parameters (per_day)"
    legend_constants[42] = "p28 in component model_parameters (per_day)"
    legend_constants[43] = "pAA17 in component model_parameters (dimensionless)"
    legend_states[18] = "c018 in component c018 (micromole_per_litre)"
    legend_constants[44] = "pAA18 in component model_parameters (dimensionless)"
    legend_states[19] = "c019 in component c019 (micromole_per_litre)"
    legend_algebraic[2] = "gt in component gt (dimensionless)"
    legend_constants[45] = "p11 in component model_parameters (per_day)"
    legend_constants[46] = "s11 in component c19 (dimensionless)"
    legend_states[20] = "c11 in component c11 (micromole_per_litre)"
    legend_constants[47] = "p12 in component model_parameters (per_day)"
    legend_constants[48] = "s12 in component c19 (dimensionless)"
    legend_states[21] = "c12 in component c12 (micromole_per_litre)"
    legend_constants[49] = "p13 in component model_parameters (per_day)"
    legend_constants[50] = "s13 in component c19 (dimensionless)"
    legend_states[22] = "c13 in component c13 (micromole_per_litre)"
    legend_constants[51] = "p10 in component model_parameters (per_day)"
    legend_constants[52] = "s10 in component c19 (dimensionless)"
    legend_states[23] = "c10 in component c10 (micromole_per_litre)"
    legend_states[24] = "c20 in component c20 (micromole_per_litre)"
    legend_states[25] = "c020 in component c020 (micromole_per_litre)"
    legend_states[26] = "c21 in component c21 (micromole_per_litre)"
    legend_states[27] = "c021 in component c021 (micromole_per_litre)"
    legend_states[28] = "c22 in component c22 (micromole_per_litre)"
    legend_states[29] = "c022 in component c022 (micromole_per_litre)"
    legend_states[30] = "c23 in component c23 (micromole_per_litre)"
    legend_states[31] = "c023 in component c023 (micromole_per_litre)"
    legend_constants[53] = "p31 in component model_parameters (per_day)"
    legend_states[32] = "c24 in component c24 (micromole_per_litre)"
    legend_states[33] = "c024 in component c024 (micromole_per_litre)"
    legend_constants[54] = "p30 in component model_parameters (dimensionless)"
    legend_constants[55] = "p32 in component model_parameters (per_day)"
    legend_constants[56] = "p24 in component model_parameters (per_day)"
    legend_constants[57] = "p33 in component pp (micromole_per_litre_per_day)"
    legend_constants[58] = "t1 in component pp (day)"
    legend_constants[59] = "p27 in component model_parameters (dimensionless)"
    legend_states[34] = "c01 in component c01 (micromole_per_litre)"
    legend_constants[60] = "pAA1 in component model_parameters (dimensionless)"
    legend_states[35] = "c02 in component c02 (micromole_per_litre)"
    legend_constants[61] = "pAA2 in component model_parameters (dimensionless)"
    legend_states[36] = "c03 in component c03 (micromole_per_litre)"
    legend_constants[62] = "pAA3 in component model_parameters (dimensionless)"
    legend_states[37] = "c04 in component c04 (micromole_per_litre)"
    legend_constants[63] = "pAA4 in component model_parameters (dimensionless)"
    legend_states[38] = "c05 in component c05 (micromole_per_litre)"
    legend_constants[64] = "pAA5 in component model_parameters (dimensionless)"
    legend_states[39] = "c06 in component c06 (micromole_per_litre)"
    legend_constants[65] = "pAA6 in component model_parameters (dimensionless)"
    legend_states[40] = "c07 in component c07 (micromole_per_litre)"
    legend_constants[66] = "pAA7 in component model_parameters (dimensionless)"
    legend_states[41] = "c08 in component c08 (micromole_per_litre)"
    legend_constants[67] = "pAA8 in component model_parameters (dimensionless)"
    legend_states[42] = "c09 in component c09 (micromole_per_litre)"
    legend_constants[68] = "pAA9 in component model_parameters (dimensionless)"
    legend_states[43] = "c010 in component c010 (micromole_per_litre)"
    legend_constants[69] = "pAA10 in component model_parameters (dimensionless)"
    legend_states[44] = "c011 in component c011 (micromole_per_litre)"
    legend_constants[70] = "pAA11 in component model_parameters (dimensionless)"
    legend_states[45] = "c012 in component c012 (micromole_per_litre)"
    legend_constants[71] = "pAA12 in component model_parameters (dimensionless)"
    legend_states[46] = "c013 in component c013 (micromole_per_litre)"
    legend_constants[72] = "pAA13 in component model_parameters (dimensionless)"
    legend_states[47] = "c014 in component c014 (micromole_per_litre)"
    legend_constants[73] = "pAA14 in component model_parameters (dimensionless)"
    legend_constants[74] = "V1 in component model_parameters (litre)"
    legend_constants[75] = "CA1 in component c01 (micromole_per_litre)"
    legend_constants[76] = "CB1 in component model_parameters (micromole_per_litre)"
    legend_constants[77] = "CA2 in component c02 (micromole_per_litre)"
    legend_constants[78] = "CA3 in component c03 (micromole_per_litre)"
    legend_constants[79] = "CA4 in component c04 (micromole_per_litre)"
    legend_constants[80] = "CA5 in component c05 (micromole_per_litre)"
    legend_constants[81] = "CA6 in component c06 (micromole_per_litre)"
    legend_constants[82] = "CA7 in component c07 (micromole_per_litre)"
    legend_constants[83] = "CA8 in component c08 (micromole_per_litre)"
    legend_constants[84] = "CA9 in component c09 (micromole_per_litre)"
    legend_constants[85] = "CA10 in component c010 (micromole_per_litre)"
    legend_constants[86] = "CA11 in component c011 (micromole_per_litre)"
    legend_constants[87] = "CA12 in component c012 (micromole_per_litre)"
    legend_constants[88] = "CA13 in component c013 (micromole_per_litre)"
    legend_constants[89] = "CA14 in component c014 (micromole_per_litre)"
    legend_constants[90] = "CA15 in component c015 (micromole_per_litre)"
    legend_constants[91] = "CB15 in component model_parameters (micromole_per_litre)"
    legend_constants[92] = "CA16 in component c016 (micromole_per_litre)"
    legend_constants[93] = "CA17 in component c017 (micromole_per_litre)"
    legend_constants[94] = "CA18 in component c018 (micromole_per_litre)"
    legend_constants[95] = "CA19 in component c019 (micromole_per_litre)"
    legend_constants[96] = "CA20 in component c020 (micromole_per_litre)"
    legend_constants[97] = "CA21 in component c021 (micromole_per_litre)"
    legend_constants[98] = "CA22 in component c022 (micromole_per_litre)"
    legend_constants[99] = "CA23 in component c023 (micromole_per_litre)"
    legend_constants[100] = "CA24 in component c024 (micromole_per_litre)"
    legend_rates[0] = "d/dt c15 in component c15 (micromole_per_litre)"
    legend_rates[2] = "d/dt c16 in component c16 (micromole_per_litre)"
    legend_rates[3] = "d/dt c17 in component c17 (micromole_per_litre)"
    legend_rates[7] = "d/dt c18 in component c18 (micromole_per_litre)"
    legend_rates[4] = "d/dt c19 in component c19 (micromole_per_litre)"
    legend_rates[24] = "d/dt c20 in component c20 (micromole_per_litre)"
    legend_rates[26] = "d/dt c21 in component c21 (micromole_per_litre)"
    legend_rates[28] = "d/dt c22 in component c22 (micromole_per_litre)"
    legend_rates[30] = "d/dt c23 in component c23 (micromole_per_litre)"
    legend_rates[32] = "d/dt c24 in component c24 (micromole_per_litre)"
    legend_rates[8] = "d/dt c1 in component c1 (micromole_per_litre)"
    legend_rates[9] = "d/dt c2 in component c2 (micromole_per_litre)"
    legend_rates[10] = "d/dt c3 in component c3 (micromole_per_litre)"
    legend_rates[11] = "d/dt c4 in component c4 (micromole_per_litre)"
    legend_rates[12] = "d/dt c5 in component c5 (micromole_per_litre)"
    legend_rates[13] = "d/dt c6 in component c6 (micromole_per_litre)"
    legend_rates[14] = "d/dt c7 in component c7 (micromole_per_litre)"
    legend_rates[15] = "d/dt c8 in component c8 (micromole_per_litre)"
    legend_rates[16] = "d/dt c9 in component c9 (micromole_per_litre)"
    legend_rates[23] = "d/dt c10 in component c10 (micromole_per_litre)"
    legend_rates[20] = "d/dt c11 in component c11 (micromole_per_litre)"
    legend_rates[21] = "d/dt c12 in component c12 (micromole_per_litre)"
    legend_rates[22] = "d/dt c13 in component c13 (micromole_per_litre)"
    legend_rates[17] = "d/dt c14 in component c14 (micromole_per_litre)"
    legend_rates[34] = "d/dt c01 in component c01 (micromole_per_litre)"
    legend_rates[35] = "d/dt c02 in component c02 (micromole_per_litre)"
    legend_rates[36] = "d/dt c03 in component c03 (micromole_per_litre)"
    legend_rates[37] = "d/dt c04 in component c04 (micromole_per_litre)"
    legend_rates[38] = "d/dt c05 in component c05 (micromole_per_litre)"
    legend_rates[39] = "d/dt c06 in component c06 (micromole_per_litre)"
    legend_rates[40] = "d/dt c07 in component c07 (micromole_per_litre)"
    legend_rates[41] = "d/dt c08 in component c08 (micromole_per_litre)"
    legend_rates[42] = "d/dt c09 in component c09 (micromole_per_litre)"
    legend_rates[43] = "d/dt c010 in component c010 (micromole_per_litre)"
    legend_rates[44] = "d/dt c011 in component c011 (micromole_per_litre)"
    legend_rates[45] = "d/dt c012 in component c012 (micromole_per_litre)"
    legend_rates[46] = "d/dt c013 in component c013 (micromole_per_litre)"
    legend_rates[47] = "d/dt c014 in component c014 (micromole_per_litre)"
    legend_rates[1] = "d/dt c015 in component c015 (micromole_per_litre)"
    legend_rates[5] = "d/dt c016 in component c016 (micromole_per_litre)"
    legend_rates[6] = "d/dt c017 in component c017 (micromole_per_litre)"
    legend_rates[18] = "d/dt c018 in component c018 (micromole_per_litre)"
    legend_rates[19] = "d/dt c019 in component c019 (micromole_per_litre)"
    legend_rates[25] = "d/dt c020 in component c020 (micromole_per_litre)"
    legend_rates[27] = "d/dt c021 in component c021 (micromole_per_litre)"
    legend_rates[29] = "d/dt c022 in component c022 (micromole_per_litre)"
    legend_rates[31] = "d/dt c023 in component c023 (micromole_per_litre)"
    legend_rates[33] = "d/dt c024 in component c024 (micromole_per_litre)"
    return (legend_states, legend_algebraic, legend_voi, legend_constants)

def initConsts():
    constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
    constants[0] = 3.6
    constants[1] = 1
    constants[2] = 1.2
    constants[3] = 0
    constants[4] = 3
    constants[5] = 0.024
    states[0] = 284
    constants[6] = 1191.072
    constants[7] = 0.9
    states[1] = 284
    constants[8] = 3.593
    states[2] = 168
    constants[9] = 56.73
    states[3] = 265
    constants[10] = 0.0009
    states[4] = 41
    constants[11] = 50.08
    constants[12] = 225.76
    constants[13] = 5
    constants[14] = 2.8119
    constants[15] = 0.31
    states[5] = 168
    constants[16] = 0.25
    states[6] = 265
    constants[17] = 0.0205
    states[7] = 687
    constants[18] = 0.0033
    constants[19] = 40.68
    constants[20] = 1.816
    constants[21] = 1
    states[8] = 2258
    constants[22] = 16.29
    constants[23] = 3
    states[9] = 940
    constants[24] = 100
    constants[25] = 4
    states[10] = 2604
    constants[26] = 1.147
    constants[27] = 1
    states[11] = 987
    constants[28] = 15.99
    constants[29] = 2
    states[12] = 464
    constants[30] = 4.148
    constants[31] = 1
    states[13] = 1300
    constants[32] = 3.723
    constants[33] = 1
    states[14] = 500
    constants[34] = 1.573
    constants[35] = 1
    states[15] = 1825
    constants[36] = 1.135
    constants[37] = 1
    states[16] = 2215
    constants[38] = 9.326
    constants[39] = 1
    states[17] = 859
    constants[40] = 0
    constants[41] = 0
    constants[42] = 38.05
    constants[43] = 1.27
    states[18] = 687
    constants[44] = 0.06
    states[19] = 41
    constants[45] = 34.53
    constants[46] = 1
    states[20] = 259
    constants[47] = 10.8
    constants[48] = 1
    states[21] = 861
    constants[49] = 10.76
    constants[50] = 1
    states[22] = 1957
    constants[51] = 0.192
    constants[52] = 1
    states[23] = 327
    states[24] = 0
    states[25] = 0
    states[26] = 0
    states[27] = 0
    states[28] = 0
    states[29] = 0
    states[30] = 8159
    states[31] = 0
    constants[53] = 1.31
    states[32] = 8159
    states[33] = 0
    constants[54] = 0
    constants[55] = 478.7
    constants[56] = 0.5474
    constants[57] = 0
    constants[58] = 1
    constants[59] = 0.6804
    states[34] = 2258
    constants[60] = 1.4
    states[35] = 940
    constants[61] = 0.18
    states[36] = 2604
    constants[62] = 0
    states[37] = 987
    constants[63] = 1.19
    states[38] = 464
    constants[64] = 0.45
    states[39] = 1300
    constants[65] = 1.05
    states[40] = 500
    constants[66] = 0.62
    states[41] = 1825
    constants[67] = 1.15
    states[42] = 2215
    constants[68] = 1
    states[43] = 327
    constants[69] = 0.15
    states[44] = 259
    constants[70] = 0.77
    states[45] = 861
    constants[71] = 1.55
    states[46] = 1957
    constants[72] = 0.78
    states[47] = 859
    constants[73] = 1.81
    constants[74] = 0.6
    constants[75] = 2258
    constants[76] = 0
    constants[77] = 940
    constants[78] = 2604
    constants[79] = 987
    constants[80] = 464
    constants[81] = 1300
    constants[82] = 500
    constants[83] = 1825
    constants[84] = 2215
    constants[85] = 327
    constants[86] = 259
    constants[87] = 861
    constants[88] = 1957
    constants[89] = 859
    constants[90] = 284
    constants[91] = 1500
    constants[92] = 168
    constants[93] = 265
    constants[94] = 687
    constants[95] = 41
    constants[96] = 0
    constants[97] = 0
    constants[98] = 0
    constants[99] = 0
    constants[100] = 0
    return (states, constants)

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

def computeAlgebraic(constants, states, voi):
    algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
    states = array(states)
    voi = array(voi)
    algebraic[1] = custom_piecewise([less(voi , constants[58]), constants[57] , True, 0.00000])
    algebraic[2] = custom_piecewise([less(voi , 3.00000), 0.00000 , True, constants[59]])
    algebraic[0] = custom_piecewise([greater_equal(voi , constants[1]), constants[2] , True, constants[0]])
    algebraic[3] = custom_piecewise([greater_equal(voi , constants[4]), constants[5] , True, constants[3]])
    algebraic[4] = algebraic[0]+algebraic[3]
    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)