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)