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 = 68
sizeStates = 27
sizeConstants = 67
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 (second)"
legend_states[0] = "q_Ca_i in component environment (fmol)"
legend_states[1] = "q_Ca_SR in component environment (fmol)"
legend_states[2] = "q_IP3 in component environment (fmol)"
legend_states[3] = "q_Sa_000_IPR in component environment (fmol)"
legend_states[4] = "q_Sa_001_IPR in component environment (fmol)"
legend_states[5] = "q_Sa_010_IPR in component environment (fmol)"
legend_states[6] = "q_Sa_011_IPR in component environment (fmol)"
legend_states[7] = "q_Sa_100_IPR in component environment (fmol)"
legend_states[8] = "q_Sa_101_IPR in component environment (fmol)"
legend_states[9] = "q_Sa_110_IPR in component environment (fmol)"
legend_states[10] = "q_Sa_111_IPR in component environment (fmol)"
legend_states[11] = "q_Sb_000_IPR in component environment (fmol)"
legend_states[12] = "q_Sb_001_IPR in component environment (fmol)"
legend_states[13] = "q_Sb_010_IPR in component environment (fmol)"
legend_states[14] = "q_Sb_011_IPR in component environment (fmol)"
legend_states[15] = "q_Sb_100_IPR in component environment (fmol)"
legend_states[16] = "q_Sb_101_IPR in component environment (fmol)"
legend_states[17] = "q_Sb_110_IPR in component environment (fmol)"
legend_states[18] = "q_Sb_111_IPR in component environment (fmol)"
legend_states[19] = "q_Sc_000_IPR in component environment (fmol)"
legend_states[20] = "q_Sc_001_IPR in component environment (fmol)"
legend_states[21] = "q_Sc_010_IPR in component environment (fmol)"
legend_states[22] = "q_Sc_011_IPR in component environment (fmol)"
legend_states[23] = "q_Sc_100_IPR in component environment (fmol)"
legend_states[24] = "q_Sc_101_IPR in component environment (fmol)"
legend_states[25] = "q_Sc_110_IPR in component environment (fmol)"
legend_states[26] = "q_Sc_111_IPR in component environment (fmol)"
legend_algebraic[1] = "IP3_stim in component environment (fmol_per_sec)"
legend_algebraic[31] = "v_R_IPR_main in component IPR (fmol_per_sec)"
legend_algebraic[32] = "v_R1_a_IPR in component IPR (fmol_per_sec)"
legend_algebraic[33] = "v_R2_a_IPR in component IPR (fmol_per_sec)"
legend_algebraic[34] = "v_R3_a_IPR in component IPR (fmol_per_sec)"
legend_algebraic[35] = "v_R4_a_IPR in component IPR (fmol_per_sec)"
legend_algebraic[36] = "v_R5_a_IPR in component IPR (fmol_per_sec)"
legend_algebraic[37] = "v_R6_a_IPR in component IPR (fmol_per_sec)"
legend_algebraic[38] = "v_R7_a_IPR in component IPR (fmol_per_sec)"
legend_algebraic[39] = "v_R8_a_IPR in component IPR (fmol_per_sec)"
legend_algebraic[40] = "v_R9_a_IPR in component IPR (fmol_per_sec)"
legend_algebraic[41] = "v_R10_a_IPR in component IPR (fmol_per_sec)"
legend_algebraic[42] = "v_R11_a_IPR in component IPR (fmol_per_sec)"
legend_algebraic[43] = "v_R12_a_IPR in component IPR (fmol_per_sec)"
legend_algebraic[44] = "v_R1_b_IPR in component IPR (fmol_per_sec)"
legend_algebraic[45] = "v_R2_b_IPR in component IPR (fmol_per_sec)"
legend_algebraic[46] = "v_R3_b_IPR in component IPR (fmol_per_sec)"
legend_algebraic[47] = "v_R4_b_IPR in component IPR (fmol_per_sec)"
legend_algebraic[48] = "v_R5_b_IPR in component IPR (fmol_per_sec)"
legend_algebraic[49] = "v_R6_b_IPR in component IPR (fmol_per_sec)"
legend_algebraic[50] = "v_R7_b_IPR in component IPR (fmol_per_sec)"
legend_algebraic[51] = "v_R8_b_IPR in component IPR (fmol_per_sec)"
legend_algebraic[52] = "v_R9_b_IPR in component IPR (fmol_per_sec)"
legend_algebraic[53] = "v_R10_b_IPR in component IPR (fmol_per_sec)"
legend_algebraic[54] = "v_R11_b_IPR in component IPR (fmol_per_sec)"
legend_algebraic[55] = "v_R12_b_IPR in component IPR (fmol_per_sec)"
legend_algebraic[56] = "v_R1_c_IPR in component IPR (fmol_per_sec)"
legend_algebraic[57] = "v_R2_c_IPR in component IPR (fmol_per_sec)"
legend_algebraic[58] = "v_R3_c_IPR in component IPR (fmol_per_sec)"
legend_algebraic[59] = "v_R4_c_IPR in component IPR (fmol_per_sec)"
legend_algebraic[60] = "v_R5_c_IPR in component IPR (fmol_per_sec)"
legend_algebraic[61] = "v_R6_c_IPR in component IPR (fmol_per_sec)"
legend_algebraic[62] = "v_R7_c_IPR in component IPR (fmol_per_sec)"
legend_algebraic[63] = "v_R8_c_IPR in component IPR (fmol_per_sec)"
legend_algebraic[64] = "v_R9_c_IPR in component IPR (fmol_per_sec)"
legend_algebraic[65] = "v_R10_c_IPR in component IPR (fmol_per_sec)"
legend_algebraic[66] = "v_R11_c_IPR in component IPR (fmol_per_sec)"
legend_algebraic[67] = "v_R12_c_IPR in component IPR (fmol_per_sec)"
legend_algebraic[2] = "Ca_tot in component environment (fmol)"
legend_algebraic[0] = "Ca_gate in component environment (fmol)"
legend_constants[0] = "kappa_R_IPR_main in component IPR_parameters (fmol_per_sec)"
legend_constants[1] = "kappa_R1_a_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[2] = "kappa_R2_a_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[3] = "kappa_R3_a_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[4] = "kappa_R4_a_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[5] = "kappa_R5_a_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[6] = "kappa_R6_a_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[7] = "kappa_R7_a_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[8] = "kappa_R8_a_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[9] = "kappa_R9_a_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[10] = "kappa_R10_a_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[11] = "kappa_R11_a_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[12] = "kappa_R12_a_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[13] = "kappa_R1_b_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[14] = "kappa_R2_b_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[15] = "kappa_R3_b_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[16] = "kappa_R4_b_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[17] = "kappa_R5_b_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[18] = "kappa_R6_b_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[19] = "kappa_R7_b_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[20] = "kappa_R8_b_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[21] = "kappa_R9_b_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[22] = "kappa_R10_b_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[23] = "kappa_R11_b_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[24] = "kappa_R12_b_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[25] = "kappa_R1_c_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[26] = "kappa_R2_c_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[27] = "kappa_R3_c_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[28] = "kappa_R4_c_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[29] = "kappa_R5_c_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[30] = "kappa_R6_c_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[31] = "kappa_R7_c_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[32] = "kappa_R8_c_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[33] = "kappa_R9_c_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[34] = "kappa_R10_c_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[35] = "kappa_R11_c_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[36] = "kappa_R12_c_IPR in component IPR_parameters (fmol_per_sec)"
legend_constants[37] = "K_Ca_i in component IPR_parameters (per_fmol)"
legend_constants[38] = "K_Ca_SR in component IPR_parameters (per_fmol)"
legend_constants[39] = "K_IP3 in component IPR_parameters (per_fmol)"
legend_constants[40] = "K_Sa_000_IPR in component IPR_parameters (per_fmol)"
legend_constants[41] = "K_Sa_001_IPR in component IPR_parameters (per_fmol)"
legend_constants[42] = "K_Sa_010_IPR in component IPR_parameters (per_fmol)"
legend_constants[43] = "K_Sa_011_IPR in component IPR_parameters (per_fmol)"
legend_constants[44] = "K_Sa_100_IPR in component IPR_parameters (per_fmol)"
legend_constants[45] = "K_Sa_101_IPR in component IPR_parameters (per_fmol)"
legend_constants[46] = "K_Sa_110_IPR in component IPR_parameters (per_fmol)"
legend_constants[47] = "K_Sa_111_IPR in component IPR_parameters (per_fmol)"
legend_constants[48] = "K_Sb_000_IPR in component IPR_parameters (per_fmol)"
legend_constants[49] = "K_Sb_001_IPR in component IPR_parameters (per_fmol)"
legend_constants[50] = "K_Sb_010_IPR in component IPR_parameters (per_fmol)"
legend_constants[51] = "K_Sb_011_IPR in component IPR_parameters (per_fmol)"
legend_constants[52] = "K_Sb_100_IPR in component IPR_parameters (per_fmol)"
legend_constants[53] = "K_Sb_101_IPR in component IPR_parameters (per_fmol)"
legend_constants[54] = "K_Sb_110_IPR in component IPR_parameters (per_fmol)"
legend_constants[55] = "K_Sb_111_IPR in component IPR_parameters (per_fmol)"
legend_constants[56] = "K_Sc_000_IPR in component IPR_parameters (per_fmol)"
legend_constants[57] = "K_Sc_001_IPR in component IPR_parameters (per_fmol)"
legend_constants[58] = "K_Sc_010_IPR in component IPR_parameters (per_fmol)"
legend_constants[59] = "K_Sc_011_IPR in component IPR_parameters (per_fmol)"
legend_constants[60] = "K_Sc_100_IPR in component IPR_parameters (per_fmol)"
legend_constants[61] = "K_Sc_101_IPR in component IPR_parameters (per_fmol)"
legend_constants[62] = "K_Sc_110_IPR in component IPR_parameters (per_fmol)"
legend_constants[63] = "K_Sc_111_IPR in component IPR_parameters (per_fmol)"
legend_constants[64] = "R in component constants (J_per_K_per_mol)"
legend_constants[65] = "T in component constants (kelvin)"
legend_algebraic[3] = "mu_Ca_i in component IPR (J_per_mol)"
legend_algebraic[4] = "mu_Ca_SR in component IPR (J_per_mol)"
legend_algebraic[5] = "mu_IP3 in component IPR (J_per_mol)"
legend_algebraic[7] = "mu_Sa_000_IPR in component IPR (J_per_mol)"
legend_algebraic[8] = "mu_Sa_001_IPR in component IPR (J_per_mol)"
legend_algebraic[9] = "mu_Sa_010_IPR in component IPR (J_per_mol)"
legend_algebraic[10] = "mu_Sa_011_IPR in component IPR (J_per_mol)"
legend_algebraic[11] = "mu_Sa_100_IPR in component IPR (J_per_mol)"
legend_algebraic[12] = "mu_Sa_101_IPR in component IPR (J_per_mol)"
legend_algebraic[13] = "mu_Sa_110_IPR in component IPR (J_per_mol)"
legend_algebraic[14] = "mu_Sa_111_IPR in component IPR (J_per_mol)"
legend_algebraic[15] = "mu_Sb_000_IPR in component IPR (J_per_mol)"
legend_algebraic[16] = "mu_Sb_001_IPR in component IPR (J_per_mol)"
legend_algebraic[17] = "mu_Sb_010_IPR in component IPR (J_per_mol)"
legend_algebraic[18] = "mu_Sb_011_IPR in component IPR (J_per_mol)"
legend_algebraic[19] = "mu_Sb_100_IPR in component IPR (J_per_mol)"
legend_algebraic[20] = "mu_Sb_101_IPR in component IPR (J_per_mol)"
legend_algebraic[21] = "mu_Sb_110_IPR in component IPR (J_per_mol)"
legend_algebraic[22] = "mu_Sb_111_IPR in component IPR (J_per_mol)"
legend_algebraic[23] = "mu_Sc_000_IPR in component IPR (J_per_mol)"
legend_algebraic[24] = "mu_Sc_001_IPR in component IPR (J_per_mol)"
legend_algebraic[25] = "mu_Sc_010_IPR in component IPR (J_per_mol)"
legend_algebraic[26] = "mu_Sc_011_IPR in component IPR (J_per_mol)"
legend_algebraic[27] = "mu_Sc_100_IPR in component IPR (J_per_mol)"
legend_algebraic[28] = "mu_Sc_101_IPR in component IPR (J_per_mol)"
legend_algebraic[29] = "mu_Sc_110_IPR in component IPR (J_per_mol)"
legend_algebraic[30] = "mu_Sc_111_IPR in component IPR (J_per_mol)"
legend_algebraic[6] = "v_R_IPR_main_nogate in component IPR (fmol_per_sec)"
legend_constants[66] = "F in component constants (C_per_mol)"
legend_rates[0] = "d/dt q_Ca_i in component environment (fmol)"
legend_rates[1] = "d/dt q_Ca_SR in component environment (fmol)"
legend_rates[2] = "d/dt q_IP3 in component environment (fmol)"
legend_rates[3] = "d/dt q_Sa_000_IPR in component environment (fmol)"
legend_rates[4] = "d/dt q_Sa_001_IPR in component environment (fmol)"
legend_rates[5] = "d/dt q_Sa_010_IPR in component environment (fmol)"
legend_rates[6] = "d/dt q_Sa_011_IPR in component environment (fmol)"
legend_rates[7] = "d/dt q_Sa_100_IPR in component environment (fmol)"
legend_rates[8] = "d/dt q_Sa_101_IPR in component environment (fmol)"
legend_rates[9] = "d/dt q_Sa_110_IPR in component environment (fmol)"
legend_rates[10] = "d/dt q_Sa_111_IPR in component environment (fmol)"
legend_rates[11] = "d/dt q_Sb_000_IPR in component environment (fmol)"
legend_rates[12] = "d/dt q_Sb_001_IPR in component environment (fmol)"
legend_rates[13] = "d/dt q_Sb_010_IPR in component environment (fmol)"
legend_rates[14] = "d/dt q_Sb_011_IPR in component environment (fmol)"
legend_rates[15] = "d/dt q_Sb_100_IPR in component environment (fmol)"
legend_rates[16] = "d/dt q_Sb_101_IPR in component environment (fmol)"
legend_rates[17] = "d/dt q_Sb_110_IPR in component environment (fmol)"
legend_rates[18] = "d/dt q_Sb_111_IPR in component environment (fmol)"
legend_rates[19] = "d/dt q_Sc_000_IPR in component environment (fmol)"
legend_rates[20] = "d/dt q_Sc_001_IPR in component environment (fmol)"
legend_rates[21] = "d/dt q_Sc_010_IPR in component environment (fmol)"
legend_rates[22] = "d/dt q_Sc_011_IPR in component environment (fmol)"
legend_rates[23] = "d/dt q_Sc_100_IPR in component environment (fmol)"
legend_rates[24] = "d/dt q_Sc_101_IPR in component environment (fmol)"
legend_rates[25] = "d/dt q_Sc_110_IPR in component environment (fmol)"
legend_rates[26] = "d/dt q_Sc_111_IPR in component environment (fmol)"
return (legend_states, legend_algebraic, legend_voi, legend_constants)
def initConsts():
constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
states[0] = 0.007540192
states[1] = 0.763997897
states[2] = 5e-3
states[3] = 5e-6
states[4] = 5e-6
states[5] = 5e-6
states[6] = 5e-6
states[7] = 5e-6
states[8] = 5e-6
states[9] = 5e-6
states[10] = 5e-6
states[11] = 5e-6
states[12] = 5e-6
states[13] = 5e-6
states[14] = 5e-6
states[15] = 5e-6
states[16] = 5e-6
states[17] = 5e-6
states[18] = 5e-6
states[19] = 5e-6
states[20] = 5e-6
states[21] = 5e-6
states[22] = 5e-6
states[23] = 5e-6
states[24] = 5e-6
states[25] = 5e-6
states[26] = 5e-6
constants[0] = 1.61481e+08
constants[1] = 0.0180631
constants[2] = 5.45765
constants[3] = 0.0957773
constants[4] = 1.80631
constants[5] = 7.15849e-05
constants[6] = 0.00297939
constants[7] = 0.36184
constants[8] = 6505.35
constants[9] = 0.128757
constants[10] = 0.0426298
constants[11] = 0.68272
constants[12] = 205415
constants[13] = 0.0180631
constants[14] = 5.45765
constants[15] = 0.0957773
constants[16] = 1.80631
constants[17] = 7.15849e-05
constants[18] = 0.00297939
constants[19] = 0.36184
constants[20] = 6505.35
constants[21] = 0.128757
constants[22] = 0.0426298
constants[23] = 0.68272
constants[24] = 205415
constants[25] = 0.0180631
constants[26] = 5.45765
constants[27] = 0.0957773
constants[28] = 1.80631
constants[29] = 7.15849e-05
constants[30] = 0.00297939
constants[31] = 0.36184
constants[32] = 6505.35
constants[33] = 0.128757
constants[34] = 0.0426298
constants[35] = 0.68272
constants[36] = 205415
constants[37] = 617.368
constants[38] = 1004.58
constants[39] = 523369
constants[40] = 4.83415
constants[41] = 1.59995
constants[42] = 0.911695
constants[43] = 0.301742
constants[44] = 1219.81
constants[45] = 2930.79
constants[46] = 0.000253146
constants[47] = 552.73
constants[48] = 4.83415
constants[49] = 1.59995
constants[50] = 0.911695
constants[51] = 0.301742
constants[52] = 1219.81
constants[53] = 2930.79
constants[54] = 0.000253146
constants[55] = 552.73
constants[56] = 4.83415
constants[57] = 1.59995
constants[58] = 0.911695
constants[59] = 0.301742
constants[60] = 1219.81
constants[61] = 2930.79
constants[62] = 0.000253146
constants[63] = 552.73
constants[64] = 8.31
constants[65] = 310
constants[66] = 96485
return (states, constants)
def computeRates(voi, states, constants):
rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
algebraic[3] = constants[64]*constants[65]*log(constants[37]*states[0])
algebraic[4] = constants[64]*constants[65]*log(constants[38]*states[1])
algebraic[13] = constants[64]*constants[65]*log(constants[46]*states[9])
algebraic[21] = constants[64]*constants[65]*log(constants[54]*states[17])
algebraic[29] = constants[64]*constants[65]*log(constants[62]*states[25])
algebraic[31] = constants[0]*(exp((algebraic[4]+algebraic[13]+algebraic[21]+algebraic[29])/(constants[64]*constants[65]))-exp((algebraic[3]+algebraic[13]+algebraic[21]+algebraic[29])/(constants[64]*constants[65])))
rates[1] = -algebraic[31]
algebraic[7] = constants[64]*constants[65]*log(constants[40]*states[3])
algebraic[8] = constants[64]*constants[65]*log(constants[41]*states[4])
algebraic[32] = constants[1]*(exp((algebraic[7]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[8]/(constants[64]*constants[65])))
algebraic[10] = constants[64]*constants[65]*log(constants[43]*states[6])
algebraic[33] = constants[2]*(exp((algebraic[8]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[10]/(constants[64]*constants[65])))
algebraic[5] = constants[64]*constants[65]*log(constants[39]*states[2])
algebraic[12] = constants[64]*constants[65]*log(constants[45]*states[8])
algebraic[40] = constants[9]*(exp((algebraic[8]+algebraic[5])/(constants[64]*constants[65]))-exp(algebraic[12]/(constants[64]*constants[65])))
rates[4] = (algebraic[32]-algebraic[33])-algebraic[40]
algebraic[11] = constants[64]*constants[65]*log(constants[44]*states[7])
algebraic[36] = constants[5]*(exp((algebraic[11]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[12]/(constants[64]*constants[65])))
algebraic[14] = constants[64]*constants[65]*log(constants[47]*states[10])
algebraic[37] = constants[6]*(exp((algebraic[12]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[14]/(constants[64]*constants[65])))
rates[8] = (algebraic[36]-algebraic[37])+algebraic[40]
algebraic[9] = constants[64]*constants[65]*log(constants[42]*states[5])
algebraic[35] = constants[4]*(exp(algebraic[9]/(constants[64]*constants[65]))-exp((algebraic[7]+algebraic[3])/(constants[64]*constants[65])))
algebraic[41] = constants[10]*(exp(algebraic[11]/(constants[64]*constants[65]))-exp((algebraic[7]+algebraic[5])/(constants[64]*constants[65])))
rates[3] = -algebraic[32]+algebraic[35]+algebraic[41]
algebraic[39] = constants[8]*(exp(algebraic[13]/(constants[64]*constants[65]))-exp((algebraic[11]+algebraic[3])/(constants[64]*constants[65])))
rates[7] = (-algebraic[36]+algebraic[39])-algebraic[41]
algebraic[34] = constants[3]*(exp(algebraic[10]/(constants[64]*constants[65]))-exp((algebraic[9]+algebraic[3])/(constants[64]*constants[65])))
algebraic[42] = constants[11]*(exp((algebraic[10]+algebraic[5])/(constants[64]*constants[65]))-exp(algebraic[14]/(constants[64]*constants[65])))
rates[6] = (algebraic[33]-algebraic[34])-algebraic[42]
algebraic[38] = constants[7]*(exp(algebraic[14]/(constants[64]*constants[65]))-exp((algebraic[13]+algebraic[3])/(constants[64]*constants[65])))
rates[10] = (algebraic[37]-algebraic[38])+algebraic[42]
algebraic[43] = constants[12]*(exp(algebraic[13]/(constants[64]*constants[65]))-exp((algebraic[9]+algebraic[5])/(constants[64]*constants[65])))
rates[5] = (algebraic[34]-algebraic[35])+algebraic[43]
rates[9] = (algebraic[38]-algebraic[39])-algebraic[43]
algebraic[15] = constants[64]*constants[65]*log(constants[48]*states[11])
algebraic[16] = constants[64]*constants[65]*log(constants[49]*states[12])
algebraic[44] = constants[13]*(exp((algebraic[15]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[16]/(constants[64]*constants[65])))
algebraic[18] = constants[64]*constants[65]*log(constants[51]*states[14])
algebraic[45] = constants[14]*(exp((algebraic[16]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[18]/(constants[64]*constants[65])))
algebraic[20] = constants[64]*constants[65]*log(constants[53]*states[16])
algebraic[52] = constants[21]*(exp((algebraic[16]+algebraic[5])/(constants[64]*constants[65]))-exp(algebraic[20]/(constants[64]*constants[65])))
rates[12] = (algebraic[44]-algebraic[45])-algebraic[52]
algebraic[19] = constants[64]*constants[65]*log(constants[52]*states[15])
algebraic[48] = constants[17]*(exp((algebraic[19]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[20]/(constants[64]*constants[65])))
algebraic[22] = constants[64]*constants[65]*log(constants[55]*states[18])
algebraic[49] = constants[18]*(exp((algebraic[20]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[22]/(constants[64]*constants[65])))
rates[16] = (algebraic[48]-algebraic[49])+algebraic[52]
algebraic[17] = constants[64]*constants[65]*log(constants[50]*states[13])
algebraic[47] = constants[16]*(exp(algebraic[17]/(constants[64]*constants[65]))-exp((algebraic[15]+algebraic[3])/(constants[64]*constants[65])))
algebraic[53] = constants[22]*(exp(algebraic[19]/(constants[64]*constants[65]))-exp((algebraic[15]+algebraic[5])/(constants[64]*constants[65])))
rates[11] = -algebraic[44]+algebraic[47]+algebraic[53]
algebraic[51] = constants[20]*(exp(algebraic[21]/(constants[64]*constants[65]))-exp((algebraic[19]+algebraic[3])/(constants[64]*constants[65])))
rates[15] = (-algebraic[48]+algebraic[51])-algebraic[53]
algebraic[46] = constants[15]*(exp(algebraic[18]/(constants[64]*constants[65]))-exp((algebraic[17]+algebraic[3])/(constants[64]*constants[65])))
algebraic[54] = constants[23]*(exp((algebraic[18]+algebraic[5])/(constants[64]*constants[65]))-exp(algebraic[22]/(constants[64]*constants[65])))
rates[14] = (algebraic[45]-algebraic[46])-algebraic[54]
algebraic[50] = constants[19]*(exp(algebraic[22]/(constants[64]*constants[65]))-exp((algebraic[21]+algebraic[3])/(constants[64]*constants[65])))
rates[18] = (algebraic[49]-algebraic[50])+algebraic[54]
algebraic[55] = constants[24]*(exp(algebraic[21]/(constants[64]*constants[65]))-exp((algebraic[17]+algebraic[5])/(constants[64]*constants[65])))
rates[13] = (algebraic[46]-algebraic[47])+algebraic[55]
rates[17] = (algebraic[50]-algebraic[51])-algebraic[55]
algebraic[23] = constants[64]*constants[65]*log(constants[56]*states[19])
algebraic[24] = constants[64]*constants[65]*log(constants[57]*states[20])
algebraic[56] = constants[25]*(exp((algebraic[23]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[24]/(constants[64]*constants[65])))
algebraic[26] = constants[64]*constants[65]*log(constants[59]*states[22])
algebraic[57] = constants[26]*(exp((algebraic[24]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[26]/(constants[64]*constants[65])))
algebraic[25] = constants[64]*constants[65]*log(constants[58]*states[21])
algebraic[58] = constants[27]*(exp(algebraic[26]/(constants[64]*constants[65]))-exp((algebraic[25]+algebraic[3])/(constants[64]*constants[65])))
algebraic[59] = constants[28]*(exp(algebraic[25]/(constants[64]*constants[65]))-exp((algebraic[23]+algebraic[3])/(constants[64]*constants[65])))
algebraic[27] = constants[64]*constants[65]*log(constants[60]*states[23])
algebraic[28] = constants[64]*constants[65]*log(constants[61]*states[24])
algebraic[60] = constants[29]*(exp((algebraic[27]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[28]/(constants[64]*constants[65])))
algebraic[30] = constants[64]*constants[65]*log(constants[63]*states[26])
algebraic[61] = constants[30]*(exp((algebraic[28]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[30]/(constants[64]*constants[65])))
algebraic[62] = constants[31]*(exp(algebraic[30]/(constants[64]*constants[65]))-exp((algebraic[29]+algebraic[3])/(constants[64]*constants[65])))
algebraic[63] = constants[32]*(exp(algebraic[29]/(constants[64]*constants[65]))-exp((algebraic[27]+algebraic[3])/(constants[64]*constants[65])))
rates[0] = (((((((((((((((((algebraic[31]-algebraic[32])-algebraic[33])+algebraic[34]+algebraic[35])-algebraic[36])-algebraic[37])+algebraic[38]+algebraic[39])-algebraic[44])-algebraic[45])+algebraic[46]+algebraic[47])-algebraic[48])-algebraic[49])+algebraic[50]+algebraic[51])-algebraic[56])-algebraic[57])+algebraic[58]+algebraic[59])-algebraic[60])-algebraic[61])+algebraic[62]+algebraic[63]
algebraic[64] = constants[33]*(exp((algebraic[24]+algebraic[5])/(constants[64]*constants[65]))-exp(algebraic[28]/(constants[64]*constants[65])))
rates[20] = (algebraic[56]-algebraic[57])-algebraic[64]
rates[24] = (algebraic[60]-algebraic[61])+algebraic[64]
algebraic[65] = constants[34]*(exp(algebraic[27]/(constants[64]*constants[65]))-exp((algebraic[23]+algebraic[5])/(constants[64]*constants[65])))
rates[19] = -algebraic[56]+algebraic[59]+algebraic[65]
rates[23] = (-algebraic[60]+algebraic[63])-algebraic[65]
algebraic[66] = constants[35]*(exp((algebraic[26]+algebraic[5])/(constants[64]*constants[65]))-exp(algebraic[30]/(constants[64]*constants[65])))
rates[22] = (algebraic[57]-algebraic[58])-algebraic[66]
rates[26] = (algebraic[61]-algebraic[62])+algebraic[66]
algebraic[1] = custom_piecewise([greater(voi , 0.0600000) & less(voi , 0.0700000), 100.000 , greater(voi , 0.0900000) & less(voi , 0.100000), -100.000 , True, 0.00000])
algebraic[67] = constants[36]*(exp(algebraic[29]/(constants[64]*constants[65]))-exp((algebraic[25]+algebraic[5])/(constants[64]*constants[65])))
rates[2] = ((((((((((-algebraic[40]+algebraic[41])-algebraic[42])+algebraic[43])-algebraic[52])+algebraic[53])-algebraic[54])+algebraic[55])-algebraic[64])+algebraic[65])-algebraic[66])+algebraic[67]+algebraic[1]
rates[21] = (algebraic[58]-algebraic[59])+algebraic[67]
rates[25] = (algebraic[62]-algebraic[63])-algebraic[67]
return(rates)
def computeAlgebraic(constants, states, voi):
algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
states = array(states)
voi = array(voi)
algebraic[3] = constants[64]*constants[65]*log(constants[37]*states[0])
algebraic[4] = constants[64]*constants[65]*log(constants[38]*states[1])
algebraic[13] = constants[64]*constants[65]*log(constants[46]*states[9])
algebraic[21] = constants[64]*constants[65]*log(constants[54]*states[17])
algebraic[29] = constants[64]*constants[65]*log(constants[62]*states[25])
algebraic[31] = constants[0]*(exp((algebraic[4]+algebraic[13]+algebraic[21]+algebraic[29])/(constants[64]*constants[65]))-exp((algebraic[3]+algebraic[13]+algebraic[21]+algebraic[29])/(constants[64]*constants[65])))
algebraic[7] = constants[64]*constants[65]*log(constants[40]*states[3])
algebraic[8] = constants[64]*constants[65]*log(constants[41]*states[4])
algebraic[32] = constants[1]*(exp((algebraic[7]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[8]/(constants[64]*constants[65])))
algebraic[10] = constants[64]*constants[65]*log(constants[43]*states[6])
algebraic[33] = constants[2]*(exp((algebraic[8]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[10]/(constants[64]*constants[65])))
algebraic[5] = constants[64]*constants[65]*log(constants[39]*states[2])
algebraic[12] = constants[64]*constants[65]*log(constants[45]*states[8])
algebraic[40] = constants[9]*(exp((algebraic[8]+algebraic[5])/(constants[64]*constants[65]))-exp(algebraic[12]/(constants[64]*constants[65])))
algebraic[11] = constants[64]*constants[65]*log(constants[44]*states[7])
algebraic[36] = constants[5]*(exp((algebraic[11]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[12]/(constants[64]*constants[65])))
algebraic[14] = constants[64]*constants[65]*log(constants[47]*states[10])
algebraic[37] = constants[6]*(exp((algebraic[12]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[14]/(constants[64]*constants[65])))
algebraic[9] = constants[64]*constants[65]*log(constants[42]*states[5])
algebraic[35] = constants[4]*(exp(algebraic[9]/(constants[64]*constants[65]))-exp((algebraic[7]+algebraic[3])/(constants[64]*constants[65])))
algebraic[41] = constants[10]*(exp(algebraic[11]/(constants[64]*constants[65]))-exp((algebraic[7]+algebraic[5])/(constants[64]*constants[65])))
algebraic[39] = constants[8]*(exp(algebraic[13]/(constants[64]*constants[65]))-exp((algebraic[11]+algebraic[3])/(constants[64]*constants[65])))
algebraic[34] = constants[3]*(exp(algebraic[10]/(constants[64]*constants[65]))-exp((algebraic[9]+algebraic[3])/(constants[64]*constants[65])))
algebraic[42] = constants[11]*(exp((algebraic[10]+algebraic[5])/(constants[64]*constants[65]))-exp(algebraic[14]/(constants[64]*constants[65])))
algebraic[38] = constants[7]*(exp(algebraic[14]/(constants[64]*constants[65]))-exp((algebraic[13]+algebraic[3])/(constants[64]*constants[65])))
algebraic[43] = constants[12]*(exp(algebraic[13]/(constants[64]*constants[65]))-exp((algebraic[9]+algebraic[5])/(constants[64]*constants[65])))
algebraic[15] = constants[64]*constants[65]*log(constants[48]*states[11])
algebraic[16] = constants[64]*constants[65]*log(constants[49]*states[12])
algebraic[44] = constants[13]*(exp((algebraic[15]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[16]/(constants[64]*constants[65])))
algebraic[18] = constants[64]*constants[65]*log(constants[51]*states[14])
algebraic[45] = constants[14]*(exp((algebraic[16]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[18]/(constants[64]*constants[65])))
algebraic[20] = constants[64]*constants[65]*log(constants[53]*states[16])
algebraic[52] = constants[21]*(exp((algebraic[16]+algebraic[5])/(constants[64]*constants[65]))-exp(algebraic[20]/(constants[64]*constants[65])))
algebraic[19] = constants[64]*constants[65]*log(constants[52]*states[15])
algebraic[48] = constants[17]*(exp((algebraic[19]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[20]/(constants[64]*constants[65])))
algebraic[22] = constants[64]*constants[65]*log(constants[55]*states[18])
algebraic[49] = constants[18]*(exp((algebraic[20]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[22]/(constants[64]*constants[65])))
algebraic[17] = constants[64]*constants[65]*log(constants[50]*states[13])
algebraic[47] = constants[16]*(exp(algebraic[17]/(constants[64]*constants[65]))-exp((algebraic[15]+algebraic[3])/(constants[64]*constants[65])))
algebraic[53] = constants[22]*(exp(algebraic[19]/(constants[64]*constants[65]))-exp((algebraic[15]+algebraic[5])/(constants[64]*constants[65])))
algebraic[51] = constants[20]*(exp(algebraic[21]/(constants[64]*constants[65]))-exp((algebraic[19]+algebraic[3])/(constants[64]*constants[65])))
algebraic[46] = constants[15]*(exp(algebraic[18]/(constants[64]*constants[65]))-exp((algebraic[17]+algebraic[3])/(constants[64]*constants[65])))
algebraic[54] = constants[23]*(exp((algebraic[18]+algebraic[5])/(constants[64]*constants[65]))-exp(algebraic[22]/(constants[64]*constants[65])))
algebraic[50] = constants[19]*(exp(algebraic[22]/(constants[64]*constants[65]))-exp((algebraic[21]+algebraic[3])/(constants[64]*constants[65])))
algebraic[55] = constants[24]*(exp(algebraic[21]/(constants[64]*constants[65]))-exp((algebraic[17]+algebraic[5])/(constants[64]*constants[65])))
algebraic[23] = constants[64]*constants[65]*log(constants[56]*states[19])
algebraic[24] = constants[64]*constants[65]*log(constants[57]*states[20])
algebraic[56] = constants[25]*(exp((algebraic[23]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[24]/(constants[64]*constants[65])))
algebraic[26] = constants[64]*constants[65]*log(constants[59]*states[22])
algebraic[57] = constants[26]*(exp((algebraic[24]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[26]/(constants[64]*constants[65])))
algebraic[25] = constants[64]*constants[65]*log(constants[58]*states[21])
algebraic[58] = constants[27]*(exp(algebraic[26]/(constants[64]*constants[65]))-exp((algebraic[25]+algebraic[3])/(constants[64]*constants[65])))
algebraic[59] = constants[28]*(exp(algebraic[25]/(constants[64]*constants[65]))-exp((algebraic[23]+algebraic[3])/(constants[64]*constants[65])))
algebraic[27] = constants[64]*constants[65]*log(constants[60]*states[23])
algebraic[28] = constants[64]*constants[65]*log(constants[61]*states[24])
algebraic[60] = constants[29]*(exp((algebraic[27]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[28]/(constants[64]*constants[65])))
algebraic[30] = constants[64]*constants[65]*log(constants[63]*states[26])
algebraic[61] = constants[30]*(exp((algebraic[28]+algebraic[3])/(constants[64]*constants[65]))-exp(algebraic[30]/(constants[64]*constants[65])))
algebraic[62] = constants[31]*(exp(algebraic[30]/(constants[64]*constants[65]))-exp((algebraic[29]+algebraic[3])/(constants[64]*constants[65])))
algebraic[63] = constants[32]*(exp(algebraic[29]/(constants[64]*constants[65]))-exp((algebraic[27]+algebraic[3])/(constants[64]*constants[65])))
algebraic[64] = constants[33]*(exp((algebraic[24]+algebraic[5])/(constants[64]*constants[65]))-exp(algebraic[28]/(constants[64]*constants[65])))
algebraic[65] = constants[34]*(exp(algebraic[27]/(constants[64]*constants[65]))-exp((algebraic[23]+algebraic[5])/(constants[64]*constants[65])))
algebraic[66] = constants[35]*(exp((algebraic[26]+algebraic[5])/(constants[64]*constants[65]))-exp(algebraic[30]/(constants[64]*constants[65])))
algebraic[1] = custom_piecewise([greater(voi , 0.0600000) & less(voi , 0.0700000), 100.000 , greater(voi , 0.0900000) & less(voi , 0.100000), -100.000 , True, 0.00000])
algebraic[67] = constants[36]*(exp(algebraic[29]/(constants[64]*constants[65]))-exp((algebraic[25]+algebraic[5])/(constants[64]*constants[65])))
algebraic[0] = states[4]+states[8]+states[6]+states[10]+states[5]+states[6]+states[9]+states[10]+states[12]+states[16]+states[14]+states[18]+states[13]+states[14]+states[17]+states[18]+states[20]+states[24]+states[22]+states[26]+states[21]+states[22]+states[25]+states[26]
algebraic[2] = states[0]+states[1]+algebraic[0]
algebraic[6] = constants[0]*(exp(algebraic[4]/(constants[64]*constants[65]))-exp(algebraic[3]/(constants[64]*constants[65])))
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)