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 = 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)