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 = 53
sizeStates = 7
sizeConstants = 114
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_states[0] = "OBp in component OBp (pM)"
legend_constants[0] = "OBp_t0 in component OBp (pM)"
legend_algebraic[49] = "OBp_in in component OBp (flux)"
legend_algebraic[9] = "OBp_out in component OBp (flux)"
legend_constants[1] = "OBu_t0 in component model_parameters (pM)"
legend_constants[77] = "DifferOBpNormal in component OBp (first_order_rate_constant)"
legend_constants[2] = "D_OBp_t0 in component model_parameters (first_order_rate_constant)"
legend_constants[3] = "TGFbNormal in component OBp (pM)"
legend_constants[4] = "KD_TGF_beta_repress in component model_parameters (pM)"
legend_constants[81] = "ProlifOBpNormal in component OBp (first_order_rate_constant)"
legend_constants[5] = "frac_prolifOBp_vs_differOBu in component OBp (dimensionless)"
legend_constants[6] = "OBpNormal in component OBp (pM)"
legend_constants[7] = "OBp_sat in component OBp (pM)"
legend_constants[86] = "DifferOBu in component OBp (first_order_rate_constant)"
legend_algebraic[1] = "ProlifOBp in component OBp (first_order_rate_constant)"
legend_constants[8] = "D_OBu_t0 in component model_parameters (first_order_rate_constant)"
legend_constants[9] = "pd_OBp_t0 in component model_parameters (flux)"
legend_constants[82] = "Pi_TGFbeta_OBu_act_t0 in component TGF_beta (dimensionless)"
legend_algebraic[8] = "Pi_TGFbeta_OBp_rep in component TGF_beta (dimensionless)"
legend_algebraic[46] = "Pi_WNT in component Wnt (dimensionless)"
legend_constants[112] = "Pi_WNT_0 in component Wnt (dimensionless)"
legend_states[1] = "OBa in component OBa (pM)"
legend_algebraic[10] = "OBa_in in component OBa (flux)"
legend_algebraic[2] = "OBa_out in component OBa (flux)"
legend_constants[10] = "pd_OBa_t0 in component model_parameters (flux)"
legend_constants[11] = "D_OBa_t0 in component model_parameters (first_order_rate_constant)"
legend_constants[12] = "A_OBa_t0 in component model_parameters (first_order_rate_constant)"
legend_states[2] = "OCp in component OCp (pM)"
legend_algebraic[52] = "OCp_in in component OCp (flux)"
legend_algebraic[47] = "OCp_out in component OCp (flux)"
legend_constants[13] = "OCu_t0 in component model_parameters (pM)"
legend_algebraic[44] = "Pi_RANKL_act_OCp in component RANK_RANKL_OPG (dimensionless)"
legend_algebraic[50] = "Pi_RANKL_act_OCu in component RANK_RANKL_OPG (dimensionless)"
legend_constants[14] = "D_OCu_t0 in component model_parameters (first_order_rate_constant)"
legend_constants[15] = "pd_OCp_t0 in component model_parameters (flux)"
legend_constants[16] = "D_OCp_t0 in component model_parameters (first_order_rate_constant)"
legend_states[3] = "OCa in component OCa (pM)"
legend_algebraic[51] = "OCa_in in component OCa (flux)"
legend_algebraic[7] = "OCa_out in component OCa (flux)"
legend_constants[17] = "pd_OCa_t0 in component model_parameters (flux)"
legend_constants[18] = "A_OCa_t0 in component model_parameters (first_order_rate_constant)"
legend_algebraic[5] = "Pi_TGFbeta_OCa_act in component TGF_beta (dimensionless)"
legend_states[4] = "fbm in component fbm (dimensionless)"
legend_constants[19] = "K_form in component model_parameters (second_order_rate_constant)"
legend_constants[20] = "K_res in component model_parameters (second_order_rate_constant)"
legend_algebraic[0] = "dfbmdt in component fbm (first_order_rate_constant)"
legend_states[5] = "OCY in component OCY (pM)"
legend_algebraic[3] = "OCY_act in component OCY (pM)"
legend_constants[79] = "OCY_act_0 in component OCY (pM)"
legend_constants[21] = "fact_0 in component model_parameters (pM)"
legend_states[6] = "vm in component vm (dimensionless)"
legend_constants[22] = "XKAPPA in component model_parameters (first_order_rate_constant)"
legend_constants[23] = "vmmax in component model_parameters (dimensionless)"
legend_constants[80] = "PTH_tot in component PTH (pM)"
legend_constants[83] = "Pi_PTH_act in component PTH (dimensionless)"
legend_constants[84] = "Pi_PTH_rep in component PTH (dimensionless)"
legend_constants[24] = "Beta_PTH in component model_parameters (flux)"
legend_constants[25] = "P_PTH_d in component model_parameters (flux)"
legend_constants[26] = "Deg_PTH in component model_parameters (first_order_rate_constant)"
legend_constants[27] = "KD_PTH_act in component model_parameters (pM)"
legend_constants[28] = "KD_PTH_rep in component model_parameters (pM)"
legend_algebraic[36] = "RANK in component RANK_RANKL_OPG (pM)"
legend_algebraic[37] = "RANKL in component RANK_RANKL_OPG (pM)"
legend_algebraic[38] = "OPG in component RANK_RANKL_OPG (pM)"
legend_algebraic[39] = "RANKL_tot in component RANK_RANKL_OPG (pM)"
legend_constants[85] = "RANKL_max in component RANK_RANKL_OPG (pM)"
legend_algebraic[40] = "P_RANKL1 in component RANK_RANKL_OPG (flux)"
legend_algebraic[41] = "P_RANKL2 in component RANK_RANKL_OPG (flux)"
legend_algebraic[42] = "P_RANKL in component RANK_RANKL_OPG (flux)"
legend_constants[29] = "P_RANKL_d in component model_parameters (flux)"
legend_algebraic[34] = "Pi_NO_PTH_act_rep in component NO_PTH (dimensionless)"
legend_constants[30] = "N_RANK_OCp in component model_parameters (dimensionless)"
legend_constants[31] = "N_RANKL_OBp_max in component model_parameters (dimensionless)"
legend_constants[32] = "K_RANK_RANKL in component model_parameters (pM)"
legend_constants[33] = "Beta_OPG in component model_parameters (first_order_rate_constant)"
legend_constants[34] = "OPG_max in component model_parameters (pM)"
legend_constants[35] = "Deg_OPG in component model_parameters (first_order_rate_constant)"
legend_constants[36] = "Deg_OPG_RANKL in component model_parameters (first_order_rate_constant)"
legend_constants[37] = "K_OPG_RANKL in component model_parameters (pM)"
legend_constants[38] = "Beta_RANKL_OCY in component model_parameters (first_order_rate_constant)"
legend_constants[39] = "Beta_RANKL_OBp in component model_parameters (first_order_rate_constant)"
legend_constants[40] = "Deg_RANKL in component model_parameters (first_order_rate_constant)"
legend_constants[41] = "Deg_RANK_RANKL in component model_parameters (first_order_rate_constant)"
legend_constants[42] = "KD_RANKL_act_OCp in component model_parameters (pM)"
legend_constants[43] = "KD_RANKL_act_OCu in component model_parameters (pM)"
legend_algebraic[4] = "TGF_beta in component TGF_beta (pM)"
legend_constants[78] = "TGF_beta_t0 in component TGF_beta (pM)"
legend_constants[44] = "OCa_t0 in component TGF_beta (pM)"
legend_constants[45] = "Alpha in component model_parameters (dimensionless)"
legend_constants[46] = "KD_TGF_beta_activate in component model_parameters (pM)"
legend_algebraic[6] = "Pi_TGFbeta_OBu_act in component TGF_beta (dimensionless)"
legend_algebraic[45] = "Scl in component Scl (pM)"
legend_constants[111] = "Scl_0 in component Scl (pM)"
legend_algebraic[48] = "P_Scl_b in component Scl (flux)"
legend_constants[113] = "P_Scl_b_0 in component Scl (flux)"
legend_algebraic[33] = "A in component Scl (second_order_rate_constant)"
legend_constants[108] = "A_0 in component Scl (second_order_rate_constant)"
legend_algebraic[35] = "B in component Scl (first_order_rate_constant)"
legend_constants[109] = "B_0 in component Scl (first_order_rate_constant)"
legend_algebraic[43] = "C in component Scl (flux)"
legend_constants[110] = "C_0 in component Scl (flux)"
legend_algebraic[31] = "OCYprod in component Scl (flux)"
legend_constants[107] = "OCYprod_0 in component Scl (flux)"
legend_constants[89] = "Wnt in component Wnt (pM)"
legend_constants[47] = "Wnt_0 in component model_parameters (pM)"
legend_algebraic[29] = "Pi_eps_rep in component SED (dimensionless)"
legend_constants[106] = "Pi_eps_rep_stst in component SED (dimensionless)"
legend_constants[87] = "Beta_Scl in component Scl (first_order_rate_constant)"
legend_constants[48] = "Beta_Scl_0 in component model_parameters (first_order_rate_constant)"
legend_constants[49] = "KD_SclLRP5 in component model_parameters (pM)"
legend_constants[50] = "Deg_Scl in component model_parameters (first_order_rate_constant)"
legend_constants[51] = "Scl_max in component model_parameters (pM)"
legend_constants[52] = "KD_WntLRP5 in component model_parameters (pM)"
legend_constants[53] = "Deg_SclLRP5 in component model_parameters (first_order_rate_constant)"
legend_constants[54] = "LRP5perCell in component model_parameters (dimensionless)"
legend_algebraic[11] = "LRP5tot in component Scl (pM)"
legend_constants[88] = "LRP5tot_0 in component Scl (pM)"
legend_constants[55] = "P_Scl_d in component model_parameters (flux)"
legend_constants[90] = "PTH_tot in component NO_PTH (pM)"
legend_algebraic[30] = "NO_tot in component NO_PTH (pM)"
legend_constants[103] = "NO_eq in component NO_PTH (pM)"
legend_constants[91] = "Beta_NO in component NO_PTH (first_order_rate_constant)"
legend_constants[56] = "Beta_NO_0 in component model_parameters (first_order_rate_constant)"
legend_constants[57] = "P_NO_d in component model_parameters (flux)"
legend_constants[58] = "Deg_NO in component model_parameters (first_order_rate_constant)"
legend_algebraic[28] = "Pi_eps_act in component SED (dimensionless)"
legend_constants[102] = "Pi_eps_act_stst in component SED (dimensionless)"
legend_constants[59] = "NO_max in component model_parameters (pM)"
legend_constants[104] = "KD_rep_NO in component NO_PTH (pM)"
legend_constants[60] = "aa in component NO_PTH (dimensionless)"
legend_constants[92] = "Pi_PTH_act in component NO_PTH (dimensionless)"
legend_algebraic[32] = "Pi_NO_rep in component NO_PTH (dimensionless)"
legend_constants[61] = "lambda_s in component model_parameters (dimensionless)"
legend_constants[62] = "lambda_c in component model_parameters (dimensionless)"
legend_algebraic[26] = "SED in component SED (MPa)"
legend_algebraic[27] = "SED_bm in component SED (MPa)"
legend_constants[63] = "sig_macro_t0 in component model_parameters (MPa)"
legend_constants[64] = "de_sig_macro in component model_parameters (MPa)"
legend_constants[98] = "sig_macro in component SED (MPa)"
legend_constants[93] = "sig_1 in component SED (MPa)"
legend_constants[94] = "sig_2 in component SED (MPa)"
legend_constants[99] = "sig_3 in component SED (MPa)"
legend_constants[95] = "sig_4 in component SED (MPa)"
legend_constants[96] = "sig_5 in component SED (MPa)"
legend_constants[97] = "sig_6 in component SED (MPa)"
legend_algebraic[20] = "eps_1 in component SED (dimensionless)"
legend_algebraic[21] = "eps_2 in component SED (dimensionless)"
legend_algebraic[22] = "eps_3 in component SED (dimensionless)"
legend_algebraic[23] = "eps_4 in component SED (dimensionless)"
legend_algebraic[24] = "eps_5 in component SED (dimensionless)"
legend_algebraic[25] = "eps_6 in component SED (dimensionless)"
legend_algebraic[12] = "dens_tis in component SED (g_per_cm3)"
legend_constants[65] = "vo in component model_parameters (dimensionless)"
legend_constants[66] = "rho_m in component model_parameters (g_per_cm3)"
legend_constants[67] = "rho_o in component model_parameters (g_per_cm3)"
legend_algebraic[13] = "fvas in component SED (dimensionless)"
legend_algebraic[14] = "dens in component SED (g_per_cm3)"
legend_algebraic[15] = "dens_1 in component SED (dimensionless)"
legend_algebraic[16] = "E_mod in component SED (MPa)"
legend_constants[68] = "nu in component model_parameters (dimensionless)"
legend_algebraic[17] = "S_11 in component SED (per_MPa)"
legend_algebraic[18] = "S_12 in component SED (per_MPa)"
legend_algebraic[19] = "S_44 in component SED (per_MPa)"
legend_constants[69] = "omega in component model_parameters (dimensionless)"
legend_constants[100] = "valeq in component SED (dimensionless)"
legend_constants[70] = "taueq in component model_parameters (MPa)"
legend_constants[71] = "alphaAct in component model_parameters (dimensionless)"
legend_constants[72] = "alphaRep in component model_parameters (dimensionless)"
legend_constants[73] = "gammaAct in component model_parameters (dimensionless)"
legend_constants[74] = "gammaRep in component model_parameters (dimensionless)"
legend_constants[75] = "rhoAct in component model_parameters (dimensionless)"
legend_constants[76] = "rhoRep in component model_parameters (dimensionless)"
legend_constants[101] = "deltaAct in component SED (MPa)"
legend_constants[105] = "deltaRep in component SED (MPa)"
legend_rates[0] = "d/dt OBp in component OBp (pM)"
legend_rates[1] = "d/dt OBa in component OBa (pM)"
legend_rates[2] = "d/dt OCp in component OCp (pM)"
legend_rates[3] = "d/dt OCa in component OCa (pM)"
legend_rates[4] = "d/dt fbm in component fbm (dimensionless)"
legend_rates[5] = "d/dt OCY in component OCY (pM)"
legend_rates[6] = "d/dt vm in component vm (dimensionless)"
return (legend_states, legend_algebraic, legend_voi, legend_constants)
def initConsts():
constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
states[0] = 1.0993792851e-03
constants[0] = 1.0993792851e-03
constants[1] = 0.01
constants[2] = 0.185
constants[3] = 0.0001
constants[4] = 0.000175426051821094
constants[5] = 0.5
constants[6] = 0.001
constants[7] = 0.005
constants[8] = 0.166
constants[9] = 0
states[1] = 8.7256220470e-04
constants[10] = 0
constants[11] = 0.0212
constants[12] = 0.1908
states[2] = 1.3327554974e-03
constants[13] = 0.01
constants[14] = 0.0219
constants[15] = 0
constants[16] = 0.01958
states[3] = 1.7451244094e-05
constants[17] = 0
constants[18] = 10
states[4] = 25
constants[19] = 50
constants[20] = 2500
states[5] = 1.2093400628e-02
constants[21] = 0.000414138
states[6] = 4.1302790850e-01
constants[22] = 0.007
constants[23] = 0.516
constants[24] = 250
constants[25] = 0
constants[26] = 86
constants[27] = 0.65
constants[28] = 0.222581427709954
constants[29] = 0
constants[30] = 4160
constants[31] = 2703476
constants[32] = 10
constants[33] = 162490033.783568
constants[34] = 131.428571428571
constants[35] = 532608.695652174
constants[36] = 10.132471014805
constants[37] = 0.0151142857142857
constants[38] = 5660
constants[39] = 23600
constants[40] = 10.132471014805
constants[41] = 10.132471014805
constants[42] = 3.34
constants[43] = 16.7
constants[44] = 0.0001
constants[45] = 1
constants[46] = 0.000563278809675429
constants[47] = 170
constants[48] = 24000
constants[49] = 10
constants[50] = 1
constants[51] = 70
constants[52] = 1000
constants[53] = 50
constants[54] = 5
constants[55] = 0
constants[56] = 3440
constants[57] = 0
constants[58] = 0.0021
constants[59] = 200000000
constants[60] = 2
constants[61] = 0.450454
constants[62] = 0.900909
constants[63] = 0.50365
constants[64] = 0
constants[65] = 0.428571428571428
constants[66] = 3.2
constants[67] = 1.41
constants[68] = 0.3
constants[69] = 0.95
constants[70] = 0.006652
constants[71] = 1
constants[72] = 1
constants[73] = 7
constants[74] = 8.01559
constants[75] = 0
constants[76] = 0
constants[77] = (constants[2]*1.00000)/(1.00000+constants[3]/constants[4])
constants[78] = constants[45]*constants[44]
constants[79] = constants[21]*20.0000
constants[80] = (constants[24]+constants[25])/constants[26]
constants[81] = (constants[5]*constants[77])/(1.00000-constants[6]/constants[7])
constants[82] = constants[78]/(constants[46]+constants[78])
constants[83] = constants[80]/(constants[80]+constants[27])
constants[84] = 1.00000/(1.00000+constants[80]/constants[28])
constants[85] = constants[31]*constants[0]
constants[86] = (1.00000-constants[5])*constants[8]*constants[82]
constants[87] = constants[48]
constants[88] = constants[54]*constants[0]
constants[89] = constants[47]
constants[90] = (constants[24]+constants[25])/constants[26]
constants[91] = constants[56]
constants[92] = constants[90]/(constants[90]+constants[27])
constants[93] = 0.00000
constants[94] = 0.00000
constants[95] = 0.00000
constants[96] = 0.00000
constants[97] = 0.00000
constants[98] = constants[63]+constants[64]
constants[99] = constants[98]
constants[100] = constants[69]
constants[101] = constants[70]*(power((constants[71]-constants[100])/(constants[100]-constants[75]), 1.00000/constants[73]))
constants[102] = constants[75]+((constants[71]-constants[75])*(power(constants[70], constants[73])))/(power(constants[101], constants[73])+power(constants[70], constants[73]))
constants[103] = (constants[56]*constants[102]*constants[79])/(constants[58]+(constants[56]*constants[102]*constants[79])/constants[59])
constants[104] = constants[103]/constants[60]
constants[105] = constants[70]*(power((constants[100]-constants[76])/(constants[72]-constants[100]), 1.00000/constants[74]))
constants[106] = constants[72]-((constants[72]-constants[76])*(power(constants[70], constants[74])))/(power(constants[105], constants[74])+power(constants[70], constants[74]))
constants[107] = constants[48]*constants[79]*constants[106]
constants[108] = (1.00000/constants[49])*(constants[50]+constants[107]/constants[51])
constants[109] = (constants[108]*constants[49]*(1.00000+constants[47]/constants[52])+(constants[53]*constants[88])/constants[49])-(constants[55]+constants[107])/constants[49]
constants[110] = -(constants[55]+constants[107])*(1.00000+constants[47]/constants[52])
constants[111] = (-constants[109]+power(power(constants[109], 2.00000)-4.00000*constants[108]*constants[110], 1.0/2))/(2.00000*constants[108])
constants[112] = constants[47]/(constants[52]*(1.00000+constants[47]/constants[52]+constants[111]/constants[49]))
constants[113] = constants[107]*(1.00000-constants[111]/constants[51])
return (states, constants)
def computeRates(voi, states, constants):
rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
rates[6] = constants[22]*(constants[23]-states[6])-(states[6]*constants[19]*states[1])/states[4]
algebraic[0] = constants[19]*states[1]-constants[20]*states[3]
rates[4] = algebraic[0]
rates[5] = constants[21]*algebraic[0]
algebraic[4] = constants[45]*states[3]
algebraic[8] = 1.00000/(1.00000+algebraic[4]/constants[4])
algebraic[9] = constants[2]*algebraic[8]*states[0]
algebraic[10] = algebraic[9]+constants[10]
algebraic[2] = (constants[11]+constants[12])*states[1]
rates[1] = algebraic[10]-algebraic[2]
algebraic[1] = constants[81]*(1.00000-states[0]/constants[7])
algebraic[3] = constants[21]*states[4]
algebraic[12] = 1.00000+(constants[66]-1.00000)*states[6]+(constants[67]-1.00000)*constants[65]
algebraic[14] = (algebraic[12]*states[4])/100.000
algebraic[15] = algebraic[14]/1.00000
algebraic[16] = custom_piecewise([less(algebraic[15] , 1.20000), 2014.00*(power(algebraic[15], 2.50000)) , True, 1763.00*(power(algebraic[15], 3.20000))])
algebraic[17] = 1.00000/algebraic[16]
algebraic[18] = -constants[68]/algebraic[16]
algebraic[20] = algebraic[17]*constants[93]+algebraic[18]*constants[94]+algebraic[18]*constants[99]
algebraic[21] = algebraic[18]*constants[93]+algebraic[17]*constants[94]+algebraic[18]*constants[99]
algebraic[22] = algebraic[18]*constants[93]+algebraic[18]*constants[94]+algebraic[17]*constants[99]
algebraic[19] = (2.00000+2.00000*constants[68])/algebraic[16]
algebraic[23] = algebraic[19]*constants[95]
algebraic[24] = algebraic[19]*constants[96]
algebraic[25] = algebraic[19]*constants[97]
algebraic[26] = 0.500000*(constants[93]*algebraic[20]+constants[94]*algebraic[21]+constants[99]*algebraic[22]+constants[95]*algebraic[23]+constants[96]*algebraic[24]+constants[97]*algebraic[25])
algebraic[13] = 1.00000-states[4]/100.000
algebraic[27] = algebraic[26]/(power(1.00000-algebraic[13], 2.00000))
algebraic[29] = constants[72]-((constants[72]-constants[76])*(power(algebraic[27], constants[74])))/(power(constants[105], constants[74])+power(algebraic[27], constants[74]))
algebraic[31] = constants[87]*algebraic[3]*algebraic[29]
algebraic[33] = (1.00000/constants[49])*(constants[50]+algebraic[31]/constants[51])
algebraic[11] = constants[54]*states[0]
algebraic[35] = (algebraic[33]*constants[49]*(1.00000+constants[89]/constants[52])+(constants[53]*algebraic[11])/constants[49])-(constants[55]+algebraic[31])/constants[49]
algebraic[43] = -(constants[55]+algebraic[31])*(1.00000+constants[89]/constants[52])
algebraic[45] = (-algebraic[35]+power(power(algebraic[35], 2.00000)-4.00000*algebraic[33]*algebraic[43], 1.0/2))/(2.00000*algebraic[33])
algebraic[46] = constants[89]/(constants[52]*(1.00000+constants[89]/constants[52]+algebraic[45]/constants[49]))
algebraic[49] = constants[86]*constants[1]+(algebraic[1]*states[0]*algebraic[46])/constants[112]+constants[9]
rates[0] = algebraic[49]-algebraic[9]
algebraic[28] = constants[75]+((constants[71]-constants[75])*(power(algebraic[27], constants[73])))/(power(constants[101], constants[73])+power(algebraic[27], constants[73]))
algebraic[30] = (constants[91]*algebraic[28]*algebraic[3]+constants[57])/(constants[58]+(constants[91]*algebraic[28]*algebraic[3])/constants[59])
algebraic[32] = 1.00000/(1.00000+algebraic[30]/constants[104])
algebraic[34] = constants[61]*(constants[92]+algebraic[32])+constants[62]*constants[92]*algebraic[32]
rootfind_0(voi, constants, rates, states, algebraic)
algebraic[44] = algebraic[37]/(constants[42]+algebraic[37])
algebraic[47] = constants[16]*algebraic[44]*states[2]
algebraic[51] = algebraic[47]+constants[17]
algebraic[5] = algebraic[4]/(constants[46]+algebraic[4])
algebraic[7] = constants[18]*algebraic[5]*states[3]
rates[3] = algebraic[51]-algebraic[7]
algebraic[50] = algebraic[37]/(constants[43]+algebraic[37])
algebraic[52] = constants[14]*constants[13]*algebraic[50]+constants[15]
rates[2] = algebraic[52]-algebraic[47]
return(rates)
def computeAlgebraic(constants, states, voi):
algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
states = array(states)
voi = array(voi)
algebraic[0] = constants[19]*states[1]-constants[20]*states[3]
algebraic[4] = constants[45]*states[3]
algebraic[8] = 1.00000/(1.00000+algebraic[4]/constants[4])
algebraic[9] = constants[2]*algebraic[8]*states[0]
algebraic[10] = algebraic[9]+constants[10]
algebraic[2] = (constants[11]+constants[12])*states[1]
algebraic[1] = constants[81]*(1.00000-states[0]/constants[7])
algebraic[3] = constants[21]*states[4]
algebraic[12] = 1.00000+(constants[66]-1.00000)*states[6]+(constants[67]-1.00000)*constants[65]
algebraic[14] = (algebraic[12]*states[4])/100.000
algebraic[15] = algebraic[14]/1.00000
algebraic[16] = custom_piecewise([less(algebraic[15] , 1.20000), 2014.00*(power(algebraic[15], 2.50000)) , True, 1763.00*(power(algebraic[15], 3.20000))])
algebraic[17] = 1.00000/algebraic[16]
algebraic[18] = -constants[68]/algebraic[16]
algebraic[20] = algebraic[17]*constants[93]+algebraic[18]*constants[94]+algebraic[18]*constants[99]
algebraic[21] = algebraic[18]*constants[93]+algebraic[17]*constants[94]+algebraic[18]*constants[99]
algebraic[22] = algebraic[18]*constants[93]+algebraic[18]*constants[94]+algebraic[17]*constants[99]
algebraic[19] = (2.00000+2.00000*constants[68])/algebraic[16]
algebraic[23] = algebraic[19]*constants[95]
algebraic[24] = algebraic[19]*constants[96]
algebraic[25] = algebraic[19]*constants[97]
algebraic[26] = 0.500000*(constants[93]*algebraic[20]+constants[94]*algebraic[21]+constants[99]*algebraic[22]+constants[95]*algebraic[23]+constants[96]*algebraic[24]+constants[97]*algebraic[25])
algebraic[13] = 1.00000-states[4]/100.000
algebraic[27] = algebraic[26]/(power(1.00000-algebraic[13], 2.00000))
algebraic[29] = constants[72]-((constants[72]-constants[76])*(power(algebraic[27], constants[74])))/(power(constants[105], constants[74])+power(algebraic[27], constants[74]))
algebraic[31] = constants[87]*algebraic[3]*algebraic[29]
algebraic[33] = (1.00000/constants[49])*(constants[50]+algebraic[31]/constants[51])
algebraic[11] = constants[54]*states[0]
algebraic[35] = (algebraic[33]*constants[49]*(1.00000+constants[89]/constants[52])+(constants[53]*algebraic[11])/constants[49])-(constants[55]+algebraic[31])/constants[49]
algebraic[43] = -(constants[55]+algebraic[31])*(1.00000+constants[89]/constants[52])
algebraic[45] = (-algebraic[35]+power(power(algebraic[35], 2.00000)-4.00000*algebraic[33]*algebraic[43], 1.0/2))/(2.00000*algebraic[33])
algebraic[46] = constants[89]/(constants[52]*(1.00000+constants[89]/constants[52]+algebraic[45]/constants[49]))
algebraic[49] = constants[86]*constants[1]+(algebraic[1]*states[0]*algebraic[46])/constants[112]+constants[9]
algebraic[28] = constants[75]+((constants[71]-constants[75])*(power(algebraic[27], constants[73])))/(power(constants[101], constants[73])+power(algebraic[27], constants[73]))
algebraic[30] = (constants[91]*algebraic[28]*algebraic[3]+constants[57])/(constants[58]+(constants[91]*algebraic[28]*algebraic[3])/constants[59])
algebraic[32] = 1.00000/(1.00000+algebraic[30]/constants[104])
algebraic[34] = constants[61]*(constants[92]+algebraic[32])+constants[62]*constants[92]*algebraic[32]
algebraic[44] = algebraic[37]/(constants[42]+algebraic[37])
algebraic[47] = constants[16]*algebraic[44]*states[2]
algebraic[51] = algebraic[47]+constants[17]
algebraic[5] = algebraic[4]/(constants[46]+algebraic[4])
algebraic[7] = constants[18]*algebraic[5]*states[3]
algebraic[50] = algebraic[37]/(constants[43]+algebraic[37])
algebraic[52] = constants[14]*constants[13]*algebraic[50]+constants[15]
algebraic[6] = algebraic[4]/(constants[46]+algebraic[4])
algebraic[48] = algebraic[31]*(1.00000-algebraic[45]/constants[51])
return algebraic
initialGuess0 = None
def rootfind_0(voi, constants, rates, states, algebraic):
"""Calculate values of algebraic variables for DAE"""
from scipy.optimize import fsolve
global initialGuess0
if initialGuess0 is None: initialGuess0 = ones(7)*0.1
if not iterable(voi):
soln = fsolve(residualSN_0, initialGuess0, args=(algebraic, voi, constants, rates, states), xtol=1E-6)
initialGuess0 = soln
algebraic[36] = soln[0]
algebraic[37] = soln[1]
algebraic[38] = soln[2]
algebraic[39] = soln[3]
algebraic[40] = soln[4]
algebraic[41] = soln[5]
algebraic[42] = soln[6]
else:
for (i,t) in enumerate(voi):
soln = fsolve(residualSN_0, initialGuess0, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6)
initialGuess0 = soln
algebraic[36][i] = soln[0]
algebraic[37][i] = soln[1]
algebraic[38][i] = soln[2]
algebraic[39][i] = soln[3]
algebraic[40][i] = soln[4]
algebraic[41][i] = soln[5]
algebraic[42][i] = soln[6]
def residualSN_0(algebraicCandidate, algebraic, voi, constants, rates, states):
resid = array([0.0] * 7)
algebraic[36] = algebraicCandidate[0]
algebraic[37] = algebraicCandidate[1]
algebraic[38] = algebraicCandidate[2]
algebraic[39] = algebraicCandidate[3]
algebraic[40] = algebraicCandidate[4]
algebraic[41] = algebraicCandidate[5]
algebraic[42] = algebraicCandidate[6]
resid[0] = (algebraic[36]-(constants[30]*states[2])/(1.00000+(1.00000/constants[32])*algebraic[37]))
resid[1] = (algebraic[38]-(constants[33]*constants[84]*states[1])/((constants[33]*constants[84]*states[1])/constants[34]+constants[35]+(constants[36]/constants[37])*algebraic[37]))
resid[2] = (algebraic[39]-algebraic[37]*(1.00000+algebraic[36]/constants[32]+algebraic[38]/constants[37]))
resid[3] = (algebraic[40]-constants[38]*(1.00000-algebraic[39]/constants[85])*algebraic[3])
resid[4] = (algebraic[41]-constants[39]*algebraic[34]*(1.00000-algebraic[39]/constants[85])*states[0])
resid[5] = (algebraic[42]-(algebraic[40]+algebraic[41]))
resid[6] = (algebraic[37]-(algebraic[42]+constants[29])/(constants[40]+(constants[41]/constants[32])*algebraic[36]+(constants[36]/constants[37])*algebraic[38]))
return resid
def custom_piecewise(cases):
"""Compute result of a piecewise function"""
return select(cases[0::2],cases[1::2])
def solve_model():
"""Solve model with ODE solver"""
from scipy.integrate import ode
# Initialise constants and state variables
(init_states, constants) = initConsts()
# Set timespan to solve over
voi = linspace(0, 10, 500)
# Construct ODE object to solve
r = ode(computeRates)
r.set_integrator('vode', method='bdf', atol=1e-06, rtol=1e-06, max_step=1)
r.set_initial_value(init_states, voi[0])
r.set_f_params(constants)
# Solve model
states = array([[0.0] * len(voi)] * sizeStates)
states[:,0] = init_states
for (i,t) in enumerate(voi[1:]):
if r.successful():
r.integrate(t)
states[:,i+1] = r.y
else:
break
# Compute algebraic variables
algebraic = computeAlgebraic(constants, states, voi)
return (voi, states, algebraic)
def plot_model(voi, states, algebraic):
"""Plot variables against variable of integration"""
import pylab
(legend_states, legend_algebraic, legend_voi, legend_constants) = createLegends()
pylab.figure(1)
pylab.plot(voi,vstack((states,algebraic)).T)
pylab.xlabel(legend_voi)
pylab.legend(legend_states + legend_algebraic, loc='best')
pylab.show()
if __name__ == "__main__":
(voi, states, algebraic) = solve_model()
plot_model(voi, states, algebraic)