# Size of variable arrays: sizeAlgebraic = 71 sizeStates = 56 sizeConstants = 105 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 cell (millisecond)" legend_constants[0] = "C_m in component wal_environment (microF_per_cm2)" legend_constants[1] = "gam in component wal_environment (dimensionless)" legend_constants[2] = "R_a in component wal_environment (ohm_cm2)" legend_constants[3] = "tsi in component wal_environment (centi_metre)" legend_constants[4] = "tsi2 in component wal_environment (centi_metre)" legend_constants[5] = "tsi3 in component wal_environment (centi_metre)" legend_constants[6] = "FF in component wal_environment (C_per_mol)" legend_constants[7] = "tau_K in component wal_environment (millisecond)" legend_constants[8] = "tau_Na in component wal_environment (millisecond)" legend_constants[9] = "f_T in component wal_environment (dimensionless)" legend_constants[10] = "tau_K2 in component wal_environment (millisecond)" legend_constants[11] = "tau_Na2 in component wal_environment (millisecond)" legend_constants[12] = "I_K_rest in component wal_environment (microA_per_cm2)" legend_constants[13] = "I_Na_rest in component wal_environment (microA_per_cm2)" legend_constants[14] = "alpha_h_bar in component wal_environment (per_millisecond)" legend_constants[15] = "alpha_m_bar in component wal_environment (per_millisecond_per_millivolt)" legend_constants[16] = "alpha_n_bar in component wal_environment (per_millisecond_per_millivolt)" legend_constants[17] = "beta_h_bar in component wal_environment (per_millisecond)" legend_constants[18] = "beta_m_bar in component wal_environment (per_millisecond)" legend_constants[19] = "beta_n_bar in component wal_environment (per_millisecond)" legend_constants[20] = "V_m in component wal_environment (millivolt)" legend_constants[21] = "V_n in component wal_environment (millivolt)" legend_constants[22] = "V_h in component wal_environment (millivolt)" legend_constants[23] = "V_a in component wal_environment (millivolt)" legend_constants[24] = "V_S_inf in component wal_environment (millivolt)" legend_constants[25] = "V_h_K_inf in component wal_environment (millivolt)" legend_constants[26] = "A_a in component wal_environment (millivolt)" legend_constants[27] = "A_S_inf in component wal_environment (millivolt)" legend_constants[28] = "A_h_K_inf in component wal_environment (millivolt)" legend_constants[29] = "K_alpha_h in component wal_environment (millivolt)" legend_constants[30] = "K_beta_h in component wal_environment (millivolt)" legend_constants[31] = "K_alpha_m in component wal_environment (millivolt)" legend_constants[32] = "K_alpha_n in component wal_environment (millivolt)" legend_constants[33] = "K_beta_m in component wal_environment (millivolt)" legend_constants[34] = "K_beta_n in component wal_environment (millivolt)" legend_constants[35] = "RR in component wal_environment (milliJ_per_degreeK_per_mol)" legend_constants[36] = "TT in component wal_environment (degreeK)" legend_constants[37] = "g_Cl_bar in component wal_environment (milliS_per_cm2)" legend_constants[38] = "g_K_bar in component wal_environment (milliS_per_cm2)" legend_constants[39] = "g_Na_bar in component wal_environment (milliS_per_cm2)" legend_constants[40] = "G_K in component wal_environment (milliS_per_cm2)" legend_constants[41] = "del in component wal_environment (dimensionless)" legend_constants[42] = "K_K in component wal_environment (milliM2)" legend_constants[43] = "K_S in component wal_environment (milliM2)" legend_constants[44] = "K_m_K in component wal_environment (milliM)" legend_constants[45] = "K_m_Na in component wal_environment (milliM)" legend_constants[46] = "S_i in component wal_environment (milliM)" legend_constants[47] = "J_NaK_bar in component wal_environment (micro_mol_per_cm2_per_second)" legend_constants[48] = "V_tau in component wal_environment (millivolt)" legend_algebraic[0] = "I_T in component wal_environment (microA_per_cm2)" legend_states[0] = "vS in component wal_environment (millivolt)" legend_states[1] = "vT in component wal_environment (millivolt)" legend_algebraic[51] = "I_ionic_s in component wal_environment (microA_per_cm2)" legend_algebraic[70] = "I_ionic_t in component wal_environment (microA_per_cm2)" legend_states[2] = "K_t in component wal_environment (milliM)" legend_states[3] = "K_i in component wal_environment (milliM)" legend_states[4] = "K_e in component wal_environment (milliM)" legend_states[5] = "Na_i in component wal_environment (milliM)" legend_states[6] = "Na_t in component wal_environment (milliM)" legend_states[7] = "Na_e in component wal_environment (milliM)" legend_algebraic[13] = "E_K in component wal_environment (millivolt)" legend_algebraic[25] = "E_K_t in component wal_environment (millivolt)" legend_algebraic[26] = "Cl_i in component wal_environment (milliM)" legend_algebraic[27] = "Cl_o in component wal_environment (milliM)" legend_algebraic[28] = "Cl_i_t in component wal_environment (milliM)" legend_algebraic[29] = "Cl_o_t in component wal_environment (milliM)" legend_algebraic[30] = "J_K in component wal_environment (milliV_milliM)" legend_algebraic[31] = "J_K_t in component wal_environment (milliV_milliM)" legend_constants[49] = "eta_Cl in component wal_environment (dimensionless)" legend_constants[50] = "eta_IR in component wal_environment (dimensionless)" legend_constants[51] = "eta_DR in component wal_environment (dimensionless)" legend_constants[52] = "eta_Na in component wal_environment (dimensionless)" legend_constants[53] = "eta_NaK in component wal_environment (dimensionless)" legend_algebraic[36] = "I_Cl in component sarco_Cl_channel (microA_per_cm2)" legend_algebraic[41] = "I_IR in component sarco_IR_channel (microA_per_cm2)" legend_algebraic[43] = "I_DR in component sarco_DR_channel (microA_per_cm2)" legend_algebraic[46] = "I_Na in component sarco_Na_channel (microA_per_cm2)" legend_algebraic[50] = "I_NaK in component sarco_NaK_channel (microA_per_cm2)" legend_algebraic[55] = "I_Cl_t in component t_Cl_channel (microA_per_cm2)" legend_algebraic[60] = "I_IR_t in component t_IR_channel (microA_per_cm2)" legend_algebraic[62] = "I_DR_t in component t_DR_channel (microA_per_cm2)" legend_algebraic[65] = "I_Na_t in component t_Na_channel (microA_per_cm2)" legend_algebraic[69] = "I_NaK_t in component t_NaK_channel (microA_per_cm2)" legend_algebraic[32] = "I_HH in component wal_environment (microA_per_cm2)" legend_algebraic[33] = "a in component sarco_Cl_channel (dimensionless)" legend_algebraic[34] = "J_Cl in component sarco_Cl_channel (milliV_milliM)" legend_algebraic[35] = "g_Cl in component sarco_Cl_channel (milliS_per_cm2)" legend_algebraic[37] = "K_R in component sarco_IR_channel (milliM)" legend_algebraic[38] = "g_IR_bar in component sarco_IR_channel (milliS_per_cm2)" legend_algebraic[39] = "y in component sarco_IR_channel (dimensionless)" legend_algebraic[40] = "g_IR in component sarco_IR_channel (milliS_per_cm2)" legend_algebraic[1] = "alpha_n in component sarco_DR_channel (per_millisecond)" legend_algebraic[14] = "beta_n in component sarco_DR_channel (per_millisecond)" legend_algebraic[2] = "h_K_inf in component sarco_DR_channel (dimensionless)" legend_algebraic[15] = "tau_h_K in component sarco_DR_channel (millisecond)" legend_states[8] = "n in component sarco_DR_channel (dimensionless)" legend_states[9] = "h_K in component sarco_DR_channel (dimensionless)" legend_algebraic[42] = "g_DR in component sarco_DR_channel (milliS_per_cm2)" legend_algebraic[3] = "alpha_h in component sarco_Na_channel (per_millisecond)" legend_algebraic[16] = "beta_h in component sarco_Na_channel (per_millisecond)" legend_algebraic[4] = "alpha_m in component sarco_Na_channel (per_millisecond)" legend_algebraic[17] = "beta_m in component sarco_Na_channel (per_millisecond)" legend_algebraic[5] = "S_inf in component sarco_Na_channel (dimensionless)" legend_algebraic[18] = "tau_S in component sarco_Na_channel (millisecond)" legend_states[10] = "m in component sarco_Na_channel (dimensionless)" legend_states[11] = "h in component sarco_Na_channel (dimensionless)" legend_states[12] = "S in component sarco_Na_channel (dimensionless)" legend_algebraic[45] = "g_Na in component sarco_Na_channel (milliS_per_cm2)" legend_algebraic[44] = "J_Na in component sarco_Na_channel (milliV_milliM)" legend_algebraic[47] = "sig in component sarco_NaK_channel (dimensionless)" legend_algebraic[48] = "f1 in component sarco_NaK_channel (dimensionless)" legend_algebraic[49] = "I_NaK_bar in component sarco_NaK_channel (microA_per_cm2)" legend_algebraic[52] = "a_t in component t_Cl_channel (dimensionless)" legend_algebraic[53] = "J_Cl_t in component t_Cl_channel (milliV_milliM)" legend_algebraic[54] = "g_Cl_t in component t_Cl_channel (milliS_per_cm2)" legend_algebraic[56] = "K_R_t in component t_IR_channel (milliM)" legend_algebraic[57] = "g_IR_bar_t in component t_IR_channel (milliS_per_cm2)" legend_algebraic[58] = "y_t in component t_IR_channel (dimensionless)" legend_algebraic[59] = "g_IR_t in component t_IR_channel (milliS_per_cm2)" legend_algebraic[6] = "alpha_n_t in component t_DR_channel (per_millisecond)" legend_algebraic[19] = "beta_n_t in component t_DR_channel (per_millisecond)" legend_algebraic[7] = "h_K_inf_t in component t_DR_channel (dimensionless)" legend_algebraic[20] = "tau_h_K_t in component t_DR_channel (millisecond)" legend_states[13] = "n_t in component t_DR_channel (dimensionless)" legend_states[14] = "h_K_t in component t_DR_channel (dimensionless)" legend_algebraic[61] = "g_DR_t in component t_DR_channel (milliS_per_cm2)" legend_algebraic[8] = "alpha_h_t in component t_Na_channel (per_millisecond)" legend_algebraic[21] = "beta_h_t in component t_Na_channel (per_millisecond)" legend_algebraic[9] = "alpha_m_t in component t_Na_channel (per_millisecond)" legend_algebraic[22] = "beta_m_t in component t_Na_channel (per_millisecond)" legend_algebraic[10] = "S_inf_t in component t_Na_channel (dimensionless)" legend_algebraic[23] = "tau_S_t in component t_Na_channel (millisecond)" legend_states[15] = "m_t in component t_Na_channel (dimensionless)" legend_states[16] = "h_t in component t_Na_channel (dimensionless)" legend_states[17] = "S_t in component t_Na_channel (dimensionless)" legend_algebraic[64] = "g_Na_t in component t_Na_channel (milliS_per_cm2)" legend_algebraic[63] = "J_Na_t in component t_Na_channel (milliV_milliM)" legend_algebraic[66] = "sig_t in component t_NaK_channel (dimensionless)" legend_algebraic[67] = "f1_t in component t_NaK_channel (dimensionless)" legend_algebraic[68] = "I_NaK_bar_t in component t_NaK_channel (microA_per_cm2)" legend_states[18] = "O_0 in component sternrios (dimensionless)" legend_states[19] = "O_1 in component sternrios (dimensionless)" legend_states[20] = "O_2 in component sternrios (dimensionless)" legend_states[21] = "O_3 in component sternrios (dimensionless)" legend_states[22] = "O_4 in component sternrios (dimensionless)" legend_constants[54] = "k_L in component sternrios (per_millisecond)" legend_constants[55] = "k_Lm in component sternrios (per_millisecond)" legend_constants[56] = "f in component sternrios (dimensionless)" legend_constants[57] = "alpha1 in component sternrios (per_millisecond)" legend_constants[58] = "K in component sternrios (millivolt)" legend_constants[59] = "Vbar in component sternrios (millivolt)" legend_states[23] = "C_0 in component sternrios (dimensionless)" legend_states[24] = "C_1 in component sternrios (dimensionless)" legend_states[25] = "C_2 in component sternrios (dimensionless)" legend_states[26] = "C_3 in component sternrios (dimensionless)" legend_states[27] = "C_4 in component sternrios (dimensionless)" legend_algebraic[11] = "k_C in component sternrios (per_millisecond)" legend_algebraic[24] = "k_Cm in component sternrios (per_millisecond)" legend_constants[60] = "nu_SR in component razumova (micromolar_per_millisecond_micrometre3)" legend_constants[61] = "K_SR in component razumova (micromolar)" legend_constants[62] = "L_e in component razumova (micrometre3_per_millisecond)" legend_constants[63] = "tau_R in component razumova (micrometre3_per_millisecond)" legend_constants[64] = "tau_SR_R in component razumova (micrometre3_per_millisecond)" legend_constants[65] = "L_x in component razumova (micrometre)" legend_constants[66] = "R_R in component razumova (micrometre)" legend_constants[99] = "V_o in component razumova (micrometre3)" legend_constants[101] = "V_1 in component razumova (micrometre3)" legend_constants[102] = "V_2 in component razumova (micrometre3)" legend_constants[100] = "V_SR in component razumova (micrometre3)" legend_constants[103] = "V_SR1 in component razumova (micrometre3)" legend_constants[104] = "V_SR2 in component razumova (micrometre3)" legend_constants[67] = "k_T_on in component razumova (per_micromolar_per_millisecond)" legend_constants[68] = "k_T_off in component razumova (per_millisecond)" legend_constants[69] = "T_tot in component razumova (micromolar)" legend_constants[70] = "k_P_on in component razumova (per_micromolar_per_millisecond)" legend_constants[71] = "k_P_off in component razumova (per_millisecond)" legend_constants[72] = "P_tot in component razumova (micromolar)" legend_constants[73] = "k_Mg_on in component razumova (per_micromolar_per_millisecond)" legend_constants[74] = "k_Mg_off in component razumova (per_millisecond)" legend_constants[75] = "k_Cs_on in component razumova (per_micromolar_per_millisecond)" legend_constants[76] = "k_Cs_off in component razumova (per_millisecond)" legend_constants[77] = "Cs_tot in component razumova (micromolar)" legend_constants[78] = "k_CATP_on in component razumova (per_micromolar_per_millisecond)" legend_constants[79] = "k_CATP_off in component razumova (per_millisecond)" legend_constants[80] = "k_MATP_on in component razumova (per_micromolar_per_millisecond)" legend_constants[81] = "k_MATP_off in component razumova (per_millisecond)" legend_constants[82] = "tau_ATP in component razumova (micrometre3_per_millisecond)" legend_constants[83] = "tau_Mg in component razumova (micrometre3_per_millisecond)" legend_constants[84] = "k_0_on in component razumova (per_millisecond)" legend_constants[85] = "k_0_off in component razumova (per_millisecond)" legend_constants[86] = "k_Ca_on in component razumova (per_millisecond)" legend_constants[87] = "k_Ca_off in component razumova (per_millisecond)" legend_constants[88] = "f_o in component razumova (per_millisecond)" legend_constants[89] = "f_p in component razumova (per_millisecond)" legend_constants[90] = "h_o in component razumova (per_millisecond)" legend_constants[91] = "h_p in component razumova (per_millisecond)" legend_constants[92] = "g_o in component razumova (per_millisecond)" legend_constants[93] = "b_p in component razumova (per_millisecond)" legend_constants[94] = "k_p in component razumova (micrometre3_per_millisecond)" legend_constants[95] = "A_p in component razumova (per_milliM3_per_millisecond)" legend_constants[96] = "B_p in component razumova (per_milliM2_per_millisecond)" legend_constants[97] = "PP in component razumova (milliM2)" legend_algebraic[12] = "T_0 in component razumova (micromolar)" legend_states[28] = "Ca_1 in component razumova (micromolar)" legend_states[29] = "Ca_SR1 in component razumova (micromolar)" legend_states[30] = "Ca_2 in component razumova (micromolar)" legend_states[31] = "Ca_SR2 in component razumova (micromolar)" legend_states[32] = "Ca_T_2 in component razumova (micromolar)" legend_states[33] = "Ca_P1 in component razumova (micromolar)" legend_states[34] = "Ca_P2 in component razumova (micromolar)" legend_states[35] = "Mg_P1 in component razumova (micromolar)" legend_states[36] = "Mg_P2 in component razumova (micromolar)" legend_states[37] = "Ca_Cs1 in component razumova (micromolar)" legend_states[38] = "Ca_Cs2 in component razumova (micromolar)" legend_states[39] = "Ca_ATP1 in component razumova (micromolar)" legend_states[40] = "Ca_ATP2 in component razumova (micromolar)" legend_states[41] = "Mg_ATP1 in component razumova (micromolar)" legend_states[42] = "Mg_ATP2 in component razumova (micromolar)" legend_states[43] = "ATP1 in component razumova (micromolar)" legend_states[44] = "ATP2 in component razumova (micromolar)" legend_states[45] = "Mg1 in component razumova (micromolar)" legend_states[46] = "Mg2 in component razumova (micromolar)" legend_states[47] = "Ca_CaT2 in component razumova (micromolar)" legend_states[48] = "D_0 in component razumova (micromolar)" legend_states[49] = "D_1 in component razumova (micromolar)" legend_states[50] = "D_2 in component razumova (micromolar)" legend_states[51] = "A_1 in component razumova (micromolar)" legend_states[52] = "A_2 in component razumova (micromolar)" legend_states[53] = "P in component razumova (milliM)" legend_states[54] = "P_SR in component razumova (milliM)" legend_states[55] = "P_C_SR in component razumova (milliM)" legend_constants[98] = "i2 in component razumova (micrometre3_per_millisecond)" legend_rates[0] = "d/dt vS in component wal_environment (millivolt)" legend_rates[1] = "d/dt vT in component wal_environment (millivolt)" legend_rates[3] = "d/dt K_i in component wal_environment (milliM)" legend_rates[2] = "d/dt K_t in component wal_environment (milliM)" legend_rates[4] = "d/dt K_e in component wal_environment (milliM)" legend_rates[5] = "d/dt Na_i in component wal_environment (milliM)" legend_rates[6] = "d/dt Na_t in component wal_environment (milliM)" legend_rates[7] = "d/dt Na_e in component wal_environment (milliM)" legend_rates[8] = "d/dt n in component sarco_DR_channel (dimensionless)" legend_rates[9] = "d/dt h_K in component sarco_DR_channel (dimensionless)" legend_rates[10] = "d/dt m in component sarco_Na_channel (dimensionless)" legend_rates[11] = "d/dt h in component sarco_Na_channel (dimensionless)" legend_rates[12] = "d/dt S in component sarco_Na_channel (dimensionless)" legend_rates[13] = "d/dt n_t in component t_DR_channel (dimensionless)" legend_rates[14] = "d/dt h_K_t in component t_DR_channel (dimensionless)" legend_rates[15] = "d/dt m_t in component t_Na_channel (dimensionless)" legend_rates[16] = "d/dt h_t in component t_Na_channel (dimensionless)" legend_rates[17] = "d/dt S_t in component t_Na_channel (dimensionless)" legend_rates[23] = "d/dt C_0 in component sternrios (dimensionless)" legend_rates[18] = "d/dt O_0 in component sternrios (dimensionless)" legend_rates[24] = "d/dt C_1 in component sternrios (dimensionless)" legend_rates[19] = "d/dt O_1 in component sternrios (dimensionless)" legend_rates[25] = "d/dt C_2 in component sternrios (dimensionless)" legend_rates[20] = "d/dt O_2 in component sternrios (dimensionless)" legend_rates[26] = "d/dt C_3 in component sternrios (dimensionless)" legend_rates[21] = "d/dt O_3 in component sternrios (dimensionless)" legend_rates[27] = "d/dt C_4 in component sternrios (dimensionless)" legend_rates[22] = "d/dt O_4 in component sternrios (dimensionless)" legend_rates[28] = "d/dt Ca_1 in component razumova (micromolar)" legend_rates[29] = "d/dt Ca_SR1 in component razumova (micromolar)" legend_rates[30] = "d/dt Ca_2 in component razumova (micromolar)" legend_rates[31] = "d/dt Ca_SR2 in component razumova (micromolar)" legend_rates[32] = "d/dt Ca_T_2 in component razumova (micromolar)" legend_rates[33] = "d/dt Ca_P1 in component razumova (micromolar)" legend_rates[34] = "d/dt Ca_P2 in component razumova (micromolar)" legend_rates[35] = "d/dt Mg_P1 in component razumova (micromolar)" legend_rates[36] = "d/dt Mg_P2 in component razumova (micromolar)" legend_rates[37] = "d/dt Ca_Cs1 in component razumova (micromolar)" legend_rates[38] = "d/dt Ca_Cs2 in component razumova (micromolar)" legend_rates[39] = "d/dt Ca_ATP1 in component razumova (micromolar)" legend_rates[40] = "d/dt Ca_ATP2 in component razumova (micromolar)" legend_rates[41] = "d/dt Mg_ATP1 in component razumova (micromolar)" legend_rates[42] = "d/dt Mg_ATP2 in component razumova (micromolar)" legend_rates[43] = "d/dt ATP1 in component razumova (micromolar)" legend_rates[44] = "d/dt ATP2 in component razumova (micromolar)" legend_rates[45] = "d/dt Mg1 in component razumova (micromolar)" legend_rates[46] = "d/dt Mg2 in component razumova (micromolar)" legend_rates[47] = "d/dt Ca_CaT2 in component razumova (micromolar)" legend_rates[48] = "d/dt D_0 in component razumova (micromolar)" legend_rates[49] = "d/dt D_1 in component razumova (micromolar)" legend_rates[50] = "d/dt D_2 in component razumova (micromolar)" legend_rates[51] = "d/dt A_1 in component razumova (micromolar)" legend_rates[52] = "d/dt A_2 in component razumova (micromolar)" legend_rates[53] = "d/dt P in component razumova (milliM)" legend_rates[54] = "d/dt P_SR in component razumova (milliM)" legend_rates[55] = "d/dt P_C_SR in component razumova (milliM)" return (legend_states, legend_algebraic, legend_voi, legend_constants) def initConsts(): constants = [0.0] * sizeConstants; states = [0.0] * sizeStates; constants[0] = 0.58 constants[1] = 2.79 constants[2] = 150 constants[3] = 0.000001 constants[4] = 0.0025 constants[5] = 0.0005 constants[6] = 96485 constants[7] = 559 constants[8] = 559 constants[9] = 0.00174 constants[10] = 40229.885 constants[11] = 40229.885 constants[12] = 0.34 constants[13] = -0.43 constants[14] = 0.0081 constants[15] = 0.288 constants[16] = 0.0131 constants[17] = 4.38 constants[18] = 1.38 constants[19] = 0.067 constants[20] = -46 constants[21] = -40 constants[22] = -45 constants[23] = 70 constants[24] = -68 constants[25] = -40 constants[26] = 150 constants[27] = 7.1 constants[28] = 7.5 constants[29] = 14.7 constants[30] = 9 constants[31] = 10 constants[32] = 7 constants[33] = 18 constants[34] = 40 constants[35] = 8314.41 constants[36] = 293 constants[37] = 3.275 constants[38] = 10.8 constants[39] = 134 constants[40] = 1.85 constants[41] = 0.4 constants[42] = 950 constants[43] = 1 constants[44] = 1 constants[45] = 13 constants[46] = 10 constants[47] = 0.0001656 constants[48] = 70 states[0] = -79.974 states[1] = -80.2 states[2] = 5.9 states[3] = 150.9 states[4] = 5.9 states[5] = 12.7 states[6] = 133 states[7] = 133 constants[49] = 0.1 constants[50] = 1.0 constants[51] = 0.45 constants[52] = 0.1 constants[53] = 0.1 states[8] = 0.009466 states[9] = 0.9952 states[10] = 0.0358 states[11] = 0.4981 states[12] = 0.581 states[13] = 0.009466 states[14] = 0.9952 states[15] = 0.0358 states[16] = 0.4981 states[17] = 0.581 states[18] = 0 states[19] = 0 states[20] = 0 states[21] = 0 states[22] = 0 constants[54] = 0.002 constants[55] = 1000 constants[56] = 0.2 constants[57] = 0.2 constants[58] = 4.5 constants[59] = -20 states[23] = 1 states[24] = 0 states[25] = 0 states[26] = 0 states[27] = 0 constants[60] = 2.4375 constants[61] = 1 constants[62] = 0.00004 constants[63] = 0.75 constants[64] = 0.75 constants[65] = 1.1 constants[66] = 0.5 constants[67] = 0.0885 constants[68] = 0.115 constants[69] = 140 constants[70] = 0 constants[71] = 0 constants[72] = 1500 constants[73] = 0 constants[74] = 0 constants[75] = 0.000004 constants[76] = 0.005 constants[77] = 31000 constants[78] = 0.15 constants[79] = 30 constants[80] = 0.0015 constants[81] = 0.15 constants[82] = 0.375 constants[83] = 1.5 constants[84] = 0 constants[85] = 0.15 constants[86] = 0.15 constants[87] = 0.05 constants[88] = 0.5 constants[89] = 5 constants[90] = 0.08 constants[91] = 0.06 constants[92] = 0.04 constants[93] = 0.00000394 constants[94] = 0.00000362 constants[95] = 1 constants[96] = 0.0001 constants[97] = 6 states[28] = 0.1 states[29] = 1500 states[30] = 0.1 states[31] = 1500 states[32] = 25 states[33] = 615 states[34] = 615 states[35] = 811 states[36] = 811 states[37] = 16900 states[38] = 16900 states[39] = 0.4 states[40] = 0.4 states[41] = 7200 states[42] = 7200 states[43] = 799.6 states[44] = 799.6 states[45] = 1000 states[46] = 1000 states[47] = 3 states[48] = 0.8 states[49] = 1.2 states[50] = 3 states[51] = 0.3 states[52] = 0.23 states[53] = 0.23 states[54] = 0.23 states[55] = 0.23 constants[98] = 60 constants[99] = 0.950000*constants[65]* pi*(power(constants[66], 2.00000)) constants[100] = 0.0500000*constants[65]* pi*(power(constants[66], 2.00000)) constants[101] = 0.0100000*constants[99] constants[102] = 0.990000*constants[99] constants[103] = 0.0100000*constants[100] constants[104] = 0.990000*constants[100] return (states, constants) def computeRates(voi, states, constants): rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic rates[28] = (((((constants[98]*(states[18]+states[19]+states[20]+states[21]+states[22]))*((states[29]-states[28])/constants[101])-constants[60]*((states[28]/(states[28]+constants[61]))/constants[101]))+constants[62]*((states[29]-states[28])/constants[101]))+-constants[63]*((states[28]-states[30])/constants[101]))+-((constants[70]*states[28])*((constants[72]+-states[33])+-states[35])+-constants[71]*states[33]))+-((constants[78]*states[28])*states[43]+-constants[79]*states[39]) rates[29] = (((-(constants[98]*(states[18]+states[19]+states[20]+states[21]+states[22]))*((states[29]-states[28])/constants[103])+constants[60]*((states[28]/(states[28]+constants[61]))/constants[103]))+-constants[62]*((states[29]-states[28])/constants[103]))+-constants[64]*((states[29]-states[31])/constants[103]))+-((constants[75]*states[29])*(constants[77]-states[37])+-constants[76]*states[37]) rates[31] = (((constants[60]*((states[30]/(states[30]+constants[61]))/constants[104])+-constants[62]*((states[31]+-states[30])/constants[104]))+constants[64]*((states[29]+-states[31])/constants[104]))+-((constants[75]*states[31])*(constants[77]+-states[38])+-constants[76]*states[38]))-(1000.00/1.00000)*(constants[95]*(states[54]*(0.00100000/1.00000)*states[31]-constants[97])*(custom_piecewise([greater(states[54]*(0.00100000/1.00000)*states[31]-constants[97] , 0.00000), 1.00000 , True, 0.00000]))*(0.00100000/1.00000)*states[54]*states[31]-constants[96]*states[55]*(constants[97]-states[54]*(0.00100000/1.00000)*states[31])*(custom_piecewise([greater(constants[97]-states[54]*(0.00100000/1.00000)*states[31] , 0.00000), 1.00000 , True, 0.00000]))) rates[33] = (constants[70]*states[28])*((constants[72]+-states[33])+-states[35])+-constants[71]*states[33] rates[34] = (constants[70]*states[30])*((constants[72]+-states[34])+-states[36])+-constants[71]*states[34] rates[35] = (constants[73]*(constants[72]+-states[33]+-states[35]))*states[45]+-constants[74]*states[35] rates[36] = (constants[73]*(constants[72]+-states[34]+-states[36]))*states[46]+-constants[74]*states[36] rates[37] = (constants[75]*states[29])*(constants[77]+-states[37])+-constants[76]*states[37] rates[38] = (constants[75]*states[31])*(constants[77]+-states[38])+-constants[76]*states[38] rates[39] = ((constants[78]*states[28])*states[43]+-constants[79]*states[39])+-constants[82]*((states[39]+-states[40])/constants[101]) rates[40] = ((constants[78]*states[30])*states[44]+-constants[79]*states[40])+constants[82]*((states[39]+-states[40])/constants[102]) rates[41] = ((constants[80]*states[45])*states[43]+-constants[81]*states[41])+-constants[82]*((states[41]+-states[42])/constants[101]) rates[42] = ((constants[80]*states[46])*states[44]+-constants[81]*states[42])+constants[82]*((states[41]+-states[42])/constants[102]) rates[43] = (-((constants[78]*states[28])*states[43]+-constants[79]*states[39])+-((constants[80]*states[45])*states[43]+-constants[81]*states[41]))+-constants[82]*((states[43]+-states[44])/constants[101]) rates[44] = (-((constants[78]*states[30])*states[44]+-constants[79]*states[40])+-((constants[80]*states[46])*states[44]+-constants[81]*states[42]))+constants[82]*((states[43]+-states[44])/constants[102]) rates[45] = (-((constants[73]*(constants[72]+-states[33]+-states[35]))*states[45]+-constants[74]*states[35])+-((constants[80]*states[45])*states[43]+-constants[81]*states[41]))+-constants[83]*((states[45]+-states[46])/constants[101]) rates[46] = (-((constants[73]*(constants[72]+-states[34]+-states[36]))*states[46]+-constants[74]*states[36])+-((constants[80]*states[46])*states[44]+-constants[81]*states[42]))+constants[83]*((states[45]+-states[46])/constants[102]) rates[47] = (((constants[67]*states[30])*states[32]+-constants[68]*states[47])+-constants[86]*states[47])+constants[87]*states[50] rates[49] = ((((constants[67]*states[30]*states[48]+-constants[68]*states[49])+constants[84]*states[32])+-constants[85]*states[49])+(-constants[67]*states[30])*states[49])+constants[68]*states[50] rates[50] = (((((constants[67]*states[30]*states[49]+-constants[68]*states[50])+constants[86]*states[47])+-constants[87]*states[50])+-constants[88]*states[50])+constants[89]*states[51])+constants[92]*states[52] rates[51] = ((constants[88]*states[50]+-constants[89]*states[51])+constants[91]*states[52])+-constants[90]*states[51] rates[52] = (-constants[91]*states[52]+constants[90]*states[51])+-constants[92]*states[52] rates[53] = (0.00100000/1.00000)*(constants[90]*states[51]-constants[91]*states[52])+-1.00000*constants[93]*states[53]+-1.00000*constants[94]*((states[53]-states[54])/constants[102]) rates[54] = constants[94]*((states[53]-states[54])/constants[104])-1.00000*(constants[95]*(states[54]*(0.00100000/1.00000)*states[31]-constants[97])*(custom_piecewise([greater(states[54]*(0.00100000/1.00000)*states[31]-constants[97] , 0.00000), 1.00000 , True, 0.00000]))*(0.00100000/1.00000)*states[54]*states[31]-constants[96]*states[55]*(constants[97]-states[54]*(0.00100000/1.00000)*states[31])*(custom_piecewise([greater(constants[97]-states[54]*(0.00100000/1.00000)*states[31] , 0.00000), 1.00000 , True, 0.00000]))) rates[55] = 1.00000*(constants[95]*(states[54]*(0.00100000/1.00000)*states[31]-constants[97])*(custom_piecewise([greater(states[54]*(0.00100000/1.00000)*states[31]-constants[97] , 0.00000), 1.00000 , True, 0.00000]))*(0.00100000/1.00000)*states[54]*states[31]-constants[96]*states[55]*(constants[97]-states[54]*(0.00100000/1.00000)*states[31])*(custom_piecewise([greater(constants[97]-states[54]*(0.00100000/1.00000)*states[31] , 0.00000), 1.00000 , True, 0.00000]))) algebraic[12] = constants[69]+-states[32]+-states[47]+-states[48]+-states[49]+-states[50]+-states[51]+-states[52] rates[30] = ((((-constants[60]*((states[30]/(states[30]+constants[61]))/constants[102])+constants[62]*((states[31]+-states[30])/constants[102]))+constants[63]*((states[28]-states[30])/constants[102]))+-(((((((constants[67]*states[30]*algebraic[12]+-constants[68]*states[32])+constants[67]*states[30]*states[32])+-constants[68]*states[47])+constants[67]*states[30]*states[48])+-constants[68]*states[49])+constants[67]*states[30]*states[49])+-constants[68]*states[50]))+-((constants[70]*states[30])*(constants[72]+-states[34]+-states[36])+-constants[71]*states[34]))+-((constants[78]*states[30])*states[44]+-constants[79]*states[40]) rates[32] = (((((constants[67]*states[30])*algebraic[12]+-constants[68]*states[32])+(-constants[67]*states[30])*states[32])+constants[68]*states[47])+-constants[84]*states[32])+constants[85]*states[49] rates[48] = (((-constants[67]*states[30])*states[48]+constants[68]*states[49])+constants[84]*algebraic[12])+-constants[85]*states[48] algebraic[1] = constants[16]*((states[0]-constants[21])/(1.00000-exp(-((states[0]-constants[21])/constants[32])))) algebraic[14] = constants[19]*exp(-((states[0]-constants[21])/constants[34])) rates[8] = algebraic[1]*(1.00000-states[8])-algebraic[14]*states[8] algebraic[2] = 1.00000/(1.00000+exp((states[0]-constants[25])/constants[28])) algebraic[15] = 1000.00*exp(-((states[0]+40.0000)/25.7500)) rates[9] = (algebraic[2]-states[9])/algebraic[15] algebraic[4] = constants[15]*((states[0]-constants[20])/(1.00000-exp(-((states[0]-constants[20])/constants[31])))) algebraic[17] = constants[18]*exp(-((states[0]-constants[20])/constants[33])) rates[10] = algebraic[4]*(1.00000-states[10])-algebraic[17]*states[10] algebraic[3] = constants[14]*exp(-((states[0]-constants[22])/constants[29])) algebraic[16] = constants[17]/(1.00000+exp(-((states[0]-constants[22])/constants[30]))) rates[11] = algebraic[3]*(1.00000-states[11])-algebraic[16]*states[11] algebraic[5] = 1.00000/(1.00000+exp((states[0]-constants[24])/constants[27])) algebraic[18] = 8571.00/(0.200000+5.65000*(power((states[0]+constants[48])/100.000, 2.00000))) rates[12] = (algebraic[5]-states[12])/algebraic[18] algebraic[6] = constants[16]*((states[1]-constants[21])/(1.00000-exp(-((states[1]-constants[21])/constants[32])))) algebraic[19] = constants[19]*exp(-((states[1]-constants[21])/constants[34])) rates[13] = algebraic[6]*(1.00000-states[13])-algebraic[19]*states[13] algebraic[7] = 1.00000/(1.00000+exp((states[1]-constants[25])/constants[28])) algebraic[20] = 1.00000*exp(-((states[1]+40.0000)/25.7500)) rates[14] = (algebraic[7]-states[14])/algebraic[20] algebraic[9] = constants[15]*((states[1]-constants[20])/(1.00000-exp(-((states[1]-constants[20])/constants[31])))) algebraic[22] = constants[18]*exp(-((states[1]-constants[20])/constants[33])) rates[15] = algebraic[9]*(1.00000-states[15])-algebraic[22]*states[15] algebraic[8] = constants[14]*exp(-((states[1]-constants[22])/constants[29])) algebraic[21] = constants[17]/(1.00000+exp(-((states[1]-constants[22])/constants[30]))) rates[16] = algebraic[8]*(1.00000-states[16])-algebraic[21]*states[16] algebraic[10] = 1.00000/(1.00000+exp((states[1]-constants[24])/constants[27])) algebraic[23] = 8571.00/(0.200000+5.65000*(power((states[1]+constants[48])/100.000, 2.00000))) rates[17] = (algebraic[10]-states[17])/algebraic[23] algebraic[11] = 0.500000*constants[57]*exp((states[1]-constants[59])/(8.00000*constants[58])) algebraic[24] = 0.500000*constants[57]*exp((constants[59]-states[1])/(8.00000*constants[58])) rates[23] = -constants[54]*states[23]+constants[55]*states[18]+-4.00000*algebraic[11]*states[23]+algebraic[24]*states[24] rates[18] = constants[54]*states[23]+-constants[55]*states[18]+(-4.00000*algebraic[11]*states[18])/constants[56]+constants[56]*algebraic[24]*states[19] rates[24] = 4.00000*algebraic[11]*states[23]+-algebraic[24]*states[24]+(-constants[54]*states[24])/constants[56]+constants[56]*constants[55]*states[19]+-3.00000*algebraic[11]*states[24]+2.00000*algebraic[24]*states[25] rates[19] = (constants[54]*states[24])/constants[56]+-constants[55]*constants[56]*states[19]+(4.00000*algebraic[11]*states[18])/constants[56]+-constants[56]*algebraic[24]*states[19]+(-3.00000*algebraic[11]*states[19])/constants[56]+2.00000*constants[56]*algebraic[24]*states[20] rates[25] = 3.00000*algebraic[11]*states[24]+-2.00000*algebraic[24]*states[25]+(-constants[54]*states[25])/(power(constants[56], 2.00000))+(power(constants[56], 2.00000))*constants[55]*states[20]+-2.00000*algebraic[11]*states[25]+3.00000*algebraic[24]*states[26] rates[20] = (3.00000*algebraic[11]*states[19])/constants[56]+-2.00000*constants[56]*algebraic[24]*states[20]+(constants[54]*states[25])/(power(constants[56], 2.00000))+-constants[55]*(power(constants[56], 2.00000))*states[20]+(-2.00000*algebraic[11]*states[20])/constants[56]+3.00000*constants[56]*algebraic[24]*states[21] rates[26] = 2.00000*algebraic[11]*states[25]+-3.00000*algebraic[24]*states[26]+(-constants[54]*states[26])/(power(constants[56], 3.00000))+constants[55]*(power(constants[56], 3.00000))*states[21]+-algebraic[11]*states[26]+4.00000*algebraic[24]*states[27] rates[21] = (constants[54]*states[26])/(power(constants[56], 3.00000))+-constants[55]*(power(constants[56], 3.00000))*states[21]+(2.00000*algebraic[11]*states[20])/constants[56]+-3.00000*algebraic[24]*constants[56]*states[21]+(-algebraic[11]*states[21])/constants[56]+4.00000*constants[56]*algebraic[24]*states[22] rates[27] = algebraic[11]*states[26]+-4.00000*algebraic[24]*states[27]+(-constants[54]*states[27])/(power(constants[56], 4.00000))+constants[55]*(power(constants[56], 4.00000))*states[22] rates[22] = (algebraic[11]*states[21])/constants[56]+-4.00000*constants[56]*algebraic[24]*states[22]+(constants[54]*states[27])/(power(constants[56], 4.00000))+-constants[55]*(power(constants[56], 4.00000))*states[22] algebraic[30] = states[0]*((states[3]-states[4]*exp((-1.00000*constants[6]*states[0])/(constants[35]*constants[36])))/(1.00000-exp((-1.00000*constants[6]*states[0])/(constants[35]*constants[36])))) algebraic[13] = ((constants[35]*constants[36])/constants[6])*log(states[4]/states[3]) algebraic[37] = states[4]*exp((-constants[41]*algebraic[13])*(constants[6]/(constants[35]*constants[36]))) algebraic[38] = constants[40]*((power(algebraic[37], 2.00000))/(constants[42]+power(algebraic[37], 2.00000))) algebraic[39] = 1.00000-power(1.00000+(constants[43]*(1.00000+(power(algebraic[37], 2.00000))/constants[42]))/((power(constants[46], 2.00000))*exp((2.00000*(1.00000-constants[41])*states[0]*constants[6])/(constants[35]*constants[36]))), -1.00000) algebraic[40] = algebraic[38]*algebraic[39] algebraic[41] = algebraic[40]*(custom_piecewise([greater(algebraic[30] , 0.00000), 1.00000 , True, 0.00000]))*(algebraic[30]/50.0000) algebraic[42] = (constants[38]*(power(states[8], 4.00000)))*states[9] algebraic[43] = algebraic[42]*(algebraic[30]/50.0000) algebraic[47] = (1.00000/7.00000)*(exp(states[7]/67.3000)-1.00000) algebraic[48] = power(1.00000+0.120000*exp(-0.100000*states[0]*(constants[6]/(constants[35]*constants[36])))+0.0400000*algebraic[47]*exp(-(states[0]*(constants[6]/(constants[35]*constants[36])))), -1.00000) algebraic[49] = constants[6]*(constants[47]/((power(1.00000+constants[44]/states[4], 2.00000))*(power(1.00000+constants[45]/states[5], 3.00000)))) algebraic[50] = algebraic[49]*algebraic[48] rates[4] = (algebraic[41]+algebraic[43]+constants[12]+-2.00000*algebraic[50])/((1000.00/1.00000)*constants[6]*constants[5])+(states[2]-states[4])/constants[10] algebraic[45] = ((constants[39]*(power(states[10], 3.00000)))*states[11])*states[12] algebraic[44] = states[0]*((states[5]-states[7]*exp((-1.00000*constants[6]*states[0])/(constants[35]*constants[36])))/(1.00000-exp((-1.00000*constants[6]*states[0])/(constants[35]*constants[36])))) algebraic[46] = algebraic[45]*(algebraic[44]/75.0000) rates[7] = (algebraic[46]+constants[13]+3.00000*algebraic[50])/((1000.00/1.00000)*constants[6]*constants[5])+(states[6]-states[7])/constants[11] algebraic[0] = (1000.00/1.00000)*((states[0]-states[1])/constants[2]) algebraic[26] = 156.500/(5.00000+exp((-constants[6]*algebraic[13])/(constants[35]*constants[36]))) algebraic[27] = 156.500-5.00000*algebraic[26] algebraic[34] = states[0]*((algebraic[26]-algebraic[27]*exp((constants[6]*states[0])/(constants[35]*constants[36])))/(1.00000-exp((constants[6]*states[0])/(constants[35]*constants[36])))) algebraic[33] = 1.00000/(1.00000+exp((states[0]-constants[23])/constants[26])) algebraic[35] = constants[37]*(power(algebraic[33], 4.00000)) algebraic[36] = algebraic[35]*(algebraic[34]/45.0000) algebraic[32] = custom_piecewise([greater_equal(voi , 0.00000) & less(voi , 0.500000), 150.000 , greater_equal(voi , 50.0000) & less(voi , 50.5000), 150.000 , greater_equal(voi , 100.000) & less(voi , 100.500), 150.000 , greater_equal(voi , 150.000) & less(voi , 150.500), 150.000 , greater_equal(voi , 200.000) & less(voi , 200.500), 150.000 , greater_equal(voi , 250.000) & less(voi , 250.500), 150.000 , greater_equal(voi , 300.000) & less(voi , 300.500), 150.000 , greater_equal(voi , 350.000) & less(voi , 350.500), 150.000 , greater_equal(voi , 400.000) & less(voi , 400.500), 150.000 , True, 0.00000]) algebraic[51] = algebraic[36]+algebraic[41]+algebraic[43]+algebraic[46]+algebraic[50]+-algebraic[32] rates[0] = -((algebraic[51]+algebraic[0])/constants[0]) algebraic[31] = states[1]*((states[3]-states[2]*exp((-1.00000*constants[6]*states[1])/(constants[35]*constants[36])))/(1.00000-exp((-1.00000*constants[6]*states[1])/(constants[35]*constants[36])))) algebraic[25] = ((constants[35]*constants[36])/constants[6])*log(states[2]/states[3]) algebraic[56] = states[2]*exp((-constants[41]*algebraic[25])*(constants[6]/(constants[35]*constants[36]))) algebraic[57] = constants[40]*((power(algebraic[56], 2.00000))/(constants[42]+power(algebraic[56], 2.00000))) algebraic[58] = 1.00000-power(1.00000+(constants[43]*(1.00000+(power(algebraic[56], 2.00000))/constants[42]))/((power(constants[46], 2.00000))*exp((2.00000*(1.00000-constants[41])*states[1]*constants[6])/(constants[35]*constants[36]))), -1.00000) algebraic[59] = algebraic[57]*algebraic[58] algebraic[60] = constants[50]*algebraic[59]*(algebraic[31]/50.0000) algebraic[61] = (constants[38]*(power(states[13], 4.00000)))*states[14] algebraic[62] = constants[51]*algebraic[61]*(algebraic[31]/50.0000) algebraic[66] = (1.00000/7.00000)*(exp(states[6]/67.3000)-1.00000) algebraic[67] = power(1.00000+0.120000*exp(-0.100000*states[1]*(constants[6]/(constants[35]*constants[36])))+0.0400000*algebraic[66]*exp(-(states[1]*(constants[6]/(constants[35]*constants[36])))), -1.00000) algebraic[68] = constants[6]*(constants[47]/((power(1.00000+constants[44]/states[2], 2.00000))*(power(1.00000+constants[45]/states[5], 3.00000)))) algebraic[69] = constants[53]*algebraic[68]*algebraic[67] rates[3] = -constants[9]*((algebraic[60]+algebraic[62]+constants[12]+-2.00000*algebraic[69])/((1000.00/1.00000)*constants[6]*constants[3]))-(algebraic[41]+algebraic[43]+constants[12]+-2.00000*algebraic[50])/((1000.00/1.00000)*constants[6]*constants[4]) rates[2] = (algebraic[60]+algebraic[62]+constants[12]+-2.00000*algebraic[69])/((1000.00/1.00000)*constants[6]*constants[3])-(states[2]-states[4])/constants[7] algebraic[64] = ((constants[39]*(power(states[15], 3.00000)))*states[16])*states[17] algebraic[63] = states[1]*((states[5]-states[6]*exp((-1.00000*constants[6]*states[1])/(constants[35]*constants[36])))/(1.00000-exp((-1.00000*constants[6]*states[1])/(constants[35]*constants[36])))) algebraic[65] = constants[52]*algebraic[64]*(algebraic[63]/75.0000) rates[5] = -constants[9]*((algebraic[65]+constants[13]+3.00000*algebraic[69])/((1000.00/1.00000)*constants[6]*constants[3]))-(algebraic[46]+constants[13]+3.00000*algebraic[50])/((1000.00/1.00000)*constants[6]*constants[4]) rates[6] = (algebraic[65]+constants[13]+3.00000*algebraic[69])/((1000.00/1.00000)*constants[6]*constants[3])-(states[6]-states[7])/constants[8] algebraic[28] = 156.500/(5.00000+exp((-constants[6]*algebraic[25])/(constants[35]*constants[36]))) algebraic[29] = 156.500-5.00000*algebraic[28] algebraic[53] = states[1]*((algebraic[28]-algebraic[29]*exp((constants[6]*states[1])/(constants[35]*constants[36])))/(1.00000-exp((constants[6]*states[1])/(constants[35]*constants[36])))) algebraic[52] = 1.00000/(1.00000+exp((states[1]-constants[23])/constants[26])) algebraic[54] = constants[37]*(power(algebraic[52], 4.00000)) algebraic[55] = constants[49]*algebraic[54]*(algebraic[53]/45.0000) algebraic[70] = algebraic[55]+algebraic[60]+algebraic[62]+algebraic[65]+algebraic[69] rates[1] = -((algebraic[70]-algebraic[0]/constants[1])/constants[0]) return(rates) def computeAlgebraic(constants, states, voi): algebraic = array([[0.0] * len(voi)] * sizeAlgebraic) states = array(states) voi = array(voi) algebraic[12] = constants[69]+-states[32]+-states[47]+-states[48]+-states[49]+-states[50]+-states[51]+-states[52] algebraic[1] = constants[16]*((states[0]-constants[21])/(1.00000-exp(-((states[0]-constants[21])/constants[32])))) algebraic[14] = constants[19]*exp(-((states[0]-constants[21])/constants[34])) algebraic[2] = 1.00000/(1.00000+exp((states[0]-constants[25])/constants[28])) algebraic[15] = 1000.00*exp(-((states[0]+40.0000)/25.7500)) algebraic[4] = constants[15]*((states[0]-constants[20])/(1.00000-exp(-((states[0]-constants[20])/constants[31])))) algebraic[17] = constants[18]*exp(-((states[0]-constants[20])/constants[33])) algebraic[3] = constants[14]*exp(-((states[0]-constants[22])/constants[29])) algebraic[16] = constants[17]/(1.00000+exp(-((states[0]-constants[22])/constants[30]))) algebraic[5] = 1.00000/(1.00000+exp((states[0]-constants[24])/constants[27])) algebraic[18] = 8571.00/(0.200000+5.65000*(power((states[0]+constants[48])/100.000, 2.00000))) algebraic[6] = constants[16]*((states[1]-constants[21])/(1.00000-exp(-((states[1]-constants[21])/constants[32])))) algebraic[19] = constants[19]*exp(-((states[1]-constants[21])/constants[34])) algebraic[7] = 1.00000/(1.00000+exp((states[1]-constants[25])/constants[28])) algebraic[20] = 1.00000*exp(-((states[1]+40.0000)/25.7500)) algebraic[9] = constants[15]*((states[1]-constants[20])/(1.00000-exp(-((states[1]-constants[20])/constants[31])))) algebraic[22] = constants[18]*exp(-((states[1]-constants[20])/constants[33])) algebraic[8] = constants[14]*exp(-((states[1]-constants[22])/constants[29])) algebraic[21] = constants[17]/(1.00000+exp(-((states[1]-constants[22])/constants[30]))) algebraic[10] = 1.00000/(1.00000+exp((states[1]-constants[24])/constants[27])) algebraic[23] = 8571.00/(0.200000+5.65000*(power((states[1]+constants[48])/100.000, 2.00000))) algebraic[11] = 0.500000*constants[57]*exp((states[1]-constants[59])/(8.00000*constants[58])) algebraic[24] = 0.500000*constants[57]*exp((constants[59]-states[1])/(8.00000*constants[58])) algebraic[30] = states[0]*((states[3]-states[4]*exp((-1.00000*constants[6]*states[0])/(constants[35]*constants[36])))/(1.00000-exp((-1.00000*constants[6]*states[0])/(constants[35]*constants[36])))) algebraic[13] = ((constants[35]*constants[36])/constants[6])*log(states[4]/states[3]) algebraic[37] = states[4]*exp((-constants[41]*algebraic[13])*(constants[6]/(constants[35]*constants[36]))) algebraic[38] = constants[40]*((power(algebraic[37], 2.00000))/(constants[42]+power(algebraic[37], 2.00000))) algebraic[39] = 1.00000-power(1.00000+(constants[43]*(1.00000+(power(algebraic[37], 2.00000))/constants[42]))/((power(constants[46], 2.00000))*exp((2.00000*(1.00000-constants[41])*states[0]*constants[6])/(constants[35]*constants[36]))), -1.00000) algebraic[40] = algebraic[38]*algebraic[39] algebraic[41] = algebraic[40]*(custom_piecewise([greater(algebraic[30] , 0.00000), 1.00000 , True, 0.00000]))*(algebraic[30]/50.0000) algebraic[42] = (constants[38]*(power(states[8], 4.00000)))*states[9] algebraic[43] = algebraic[42]*(algebraic[30]/50.0000) algebraic[47] = (1.00000/7.00000)*(exp(states[7]/67.3000)-1.00000) algebraic[48] = power(1.00000+0.120000*exp(-0.100000*states[0]*(constants[6]/(constants[35]*constants[36])))+0.0400000*algebraic[47]*exp(-(states[0]*(constants[6]/(constants[35]*constants[36])))), -1.00000) algebraic[49] = constants[6]*(constants[47]/((power(1.00000+constants[44]/states[4], 2.00000))*(power(1.00000+constants[45]/states[5], 3.00000)))) algebraic[50] = algebraic[49]*algebraic[48] algebraic[45] = ((constants[39]*(power(states[10], 3.00000)))*states[11])*states[12] algebraic[44] = states[0]*((states[5]-states[7]*exp((-1.00000*constants[6]*states[0])/(constants[35]*constants[36])))/(1.00000-exp((-1.00000*constants[6]*states[0])/(constants[35]*constants[36])))) algebraic[46] = algebraic[45]*(algebraic[44]/75.0000) algebraic[0] = (1000.00/1.00000)*((states[0]-states[1])/constants[2]) algebraic[26] = 156.500/(5.00000+exp((-constants[6]*algebraic[13])/(constants[35]*constants[36]))) algebraic[27] = 156.500-5.00000*algebraic[26] algebraic[34] = states[0]*((algebraic[26]-algebraic[27]*exp((constants[6]*states[0])/(constants[35]*constants[36])))/(1.00000-exp((constants[6]*states[0])/(constants[35]*constants[36])))) algebraic[33] = 1.00000/(1.00000+exp((states[0]-constants[23])/constants[26])) algebraic[35] = constants[37]*(power(algebraic[33], 4.00000)) algebraic[36] = algebraic[35]*(algebraic[34]/45.0000) algebraic[32] = custom_piecewise([greater_equal(voi , 0.00000) & less(voi , 0.500000), 150.000 , greater_equal(voi , 50.0000) & less(voi , 50.5000), 150.000 , greater_equal(voi , 100.000) & less(voi , 100.500), 150.000 , greater_equal(voi , 150.000) & less(voi , 150.500), 150.000 , greater_equal(voi , 200.000) & less(voi , 200.500), 150.000 , greater_equal(voi , 250.000) & less(voi , 250.500), 150.000 , greater_equal(voi , 300.000) & less(voi , 300.500), 150.000 , greater_equal(voi , 350.000) & less(voi , 350.500), 150.000 , greater_equal(voi , 400.000) & less(voi , 400.500), 150.000 , True, 0.00000]) algebraic[51] = algebraic[36]+algebraic[41]+algebraic[43]+algebraic[46]+algebraic[50]+-algebraic[32] algebraic[31] = states[1]*((states[3]-states[2]*exp((-1.00000*constants[6]*states[1])/(constants[35]*constants[36])))/(1.00000-exp((-1.00000*constants[6]*states[1])/(constants[35]*constants[36])))) algebraic[25] = ((constants[35]*constants[36])/constants[6])*log(states[2]/states[3]) algebraic[56] = states[2]*exp((-constants[41]*algebraic[25])*(constants[6]/(constants[35]*constants[36]))) algebraic[57] = constants[40]*((power(algebraic[56], 2.00000))/(constants[42]+power(algebraic[56], 2.00000))) algebraic[58] = 1.00000-power(1.00000+(constants[43]*(1.00000+(power(algebraic[56], 2.00000))/constants[42]))/((power(constants[46], 2.00000))*exp((2.00000*(1.00000-constants[41])*states[1]*constants[6])/(constants[35]*constants[36]))), -1.00000) algebraic[59] = algebraic[57]*algebraic[58] algebraic[60] = constants[50]*algebraic[59]*(algebraic[31]/50.0000) algebraic[61] = (constants[38]*(power(states[13], 4.00000)))*states[14] algebraic[62] = constants[51]*algebraic[61]*(algebraic[31]/50.0000) algebraic[66] = (1.00000/7.00000)*(exp(states[6]/67.3000)-1.00000) algebraic[67] = power(1.00000+0.120000*exp(-0.100000*states[1]*(constants[6]/(constants[35]*constants[36])))+0.0400000*algebraic[66]*exp(-(states[1]*(constants[6]/(constants[35]*constants[36])))), -1.00000) algebraic[68] = constants[6]*(constants[47]/((power(1.00000+constants[44]/states[2], 2.00000))*(power(1.00000+constants[45]/states[5], 3.00000)))) algebraic[69] = constants[53]*algebraic[68]*algebraic[67] algebraic[64] = ((constants[39]*(power(states[15], 3.00000)))*states[16])*states[17] algebraic[63] = states[1]*((states[5]-states[6]*exp((-1.00000*constants[6]*states[1])/(constants[35]*constants[36])))/(1.00000-exp((-1.00000*constants[6]*states[1])/(constants[35]*constants[36])))) algebraic[65] = constants[52]*algebraic[64]*(algebraic[63]/75.0000) algebraic[28] = 156.500/(5.00000+exp((-constants[6]*algebraic[25])/(constants[35]*constants[36]))) algebraic[29] = 156.500-5.00000*algebraic[28] algebraic[53] = states[1]*((algebraic[28]-algebraic[29]*exp((constants[6]*states[1])/(constants[35]*constants[36])))/(1.00000-exp((constants[6]*states[1])/(constants[35]*constants[36])))) algebraic[52] = 1.00000/(1.00000+exp((states[1]-constants[23])/constants[26])) algebraic[54] = constants[37]*(power(algebraic[52], 4.00000)) algebraic[55] = constants[49]*algebraic[54]*(algebraic[53]/45.0000) algebraic[70] = algebraic[55]+algebraic[60]+algebraic[62]+algebraic[65]+algebraic[69] 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)