# ---------------------------------------------------------------- # This python script is the export of FCU_EC_coupling.cellml made on 4 Nov 2021 # ---------------------------------------------------------------- # Size of variable arrays: sizeAlgebraic = 199 sizeStates = 49 sizeConstants = 164 from math import * from numpy import * import matplotlib.pyplot as plt import time def createLegends(): legend_states = [""] * sizeStates legend_rates = [""] * sizeStates legend_algebraic = [""] * sizeAlgebraic legend_voi = "" legend_constants = [""] * sizeConstants legend_constants[0] = "kappa_PLBph1 in component BG_parameters (fmol_per_sec)" legend_constants[1] = "kappa_PLBph2 in component BG_parameters (fmol_per_sec)" legend_constants[2] = "kappa_PLBd1 in component BG_parameters (fmol_per_sec)" legend_constants[3] = "kappa_PLBd2 in component BG_parameters (fmol_per_sec)" legend_constants[4] = "kappa_Inh in component BG_parameters (fmol_per_sec)" legend_constants[5] = "kappa_SERCA_R1_2 in component BG_parameters (fmol_per_sec)" legend_constants[6] = "kappa_SERCA_R2_4 in component BG_parameters (fmol_per_sec)" legend_constants[7] = "kappa_SERCA_R2_2a in component BG_parameters (fmol_per_sec)" legend_constants[8] = "kappa_SERCA_R4_5 in component BG_parameters (fmol_per_sec)" legend_constants[9] = "kappa_SERCA_R5_6 in component BG_parameters (fmol_per_sec)" legend_constants[10] = "kappa_SERCA_R6_8 in component BG_parameters (fmol_per_sec)" legend_constants[11] = "kappa_SERCA_R8_9 in component BG_parameters (fmol_per_sec)" legend_constants[12] = "kappa_SERCA_R9_10 in component BG_parameters (fmol_per_sec)" legend_constants[13] = "kappa_SERCA_R10_1 in component BG_parameters (fmol_per_sec)" legend_constants[14] = "kappa_RyR in component BG_parameters (fmol_per_sec)" legend_constants[15] = "kappa_OC in component BG_parameters (fmol_per_sec)" legend_constants[16] = "kappa_CCI in component BG_parameters (fmol_per_sec)" legend_constants[17] = "kappa_CII in component BG_parameters (fmol_per_sec)" legend_constants[18] = "kappa_IO in component BG_parameters (fmol_per_sec)" legend_constants[19] = "kappa_Ca1 in component BG_parameters (fmol_per_sec)" legend_constants[20] = "kappa_Ca2 in component BG_parameters (fmol_per_sec)" legend_constants[21] = "kappa_K1 in component BG_parameters (fmol_per_sec)" legend_constants[22] = "kappa_K2 in component BG_parameters (fmol_per_sec)" legend_constants[23] = "kappa_d000 in component BG_parameters (fmol_per_sec)" legend_constants[24] = "kappa_d010 in component BG_parameters (fmol_per_sec)" legend_constants[25] = "kappa_d020 in component BG_parameters (fmol_per_sec)" legend_constants[26] = "kappa_d001 in component BG_parameters (fmol_per_sec)" legend_constants[27] = "kappa_d011 in component BG_parameters (fmol_per_sec)" legend_constants[28] = "kappa_d021 in component BG_parameters (fmol_per_sec)" legend_constants[29] = "kappa_f1_000 in component BG_parameters (fmol_per_sec)" legend_constants[30] = "kappa_f1_100 in component BG_parameters (fmol_per_sec)" legend_constants[31] = "kappa_f1_001 in component BG_parameters (fmol_per_sec)" legend_constants[32] = "kappa_f1_101 in component BG_parameters (fmol_per_sec)" legend_constants[33] = "kappa_f2_000 in component BG_parameters (fmol_per_sec)" legend_constants[34] = "kappa_f2_100 in component BG_parameters (fmol_per_sec)" legend_constants[35] = "kappa_f2_001 in component BG_parameters (fmol_per_sec)" legend_constants[36] = "kappa_f2_101 in component BG_parameters (fmol_per_sec)" legend_constants[37] = "kappa_f3_010 in component BG_parameters (fmol_per_sec)" legend_constants[38] = "kappa_f3_110 in component BG_parameters (fmol_per_sec)" legend_constants[39] = "kappa_f3_011 in component BG_parameters (fmol_per_sec)" legend_constants[40] = "kappa_f3_111 in component BG_parameters (fmol_per_sec)" legend_constants[41] = "kappa_fCa000 in component BG_parameters (fmol_per_sec)" legend_constants[42] = "kappa_fCa100 in component BG_parameters (fmol_per_sec)" legend_constants[43] = "kappa_fCa010 in component BG_parameters (fmol_per_sec)" legend_constants[44] = "kappa_fCa110 in component BG_parameters (fmol_per_sec)" legend_constants[45] = "kappa_fCa020 in component BG_parameters (fmol_per_sec)" legend_constants[46] = "kappa_fCa120 in component BG_parameters (fmol_per_sec)" legend_constants[47] = "kappa_R_TRPNCa in component BG_parameters (fmol_per_sec)" legend_constants[48] = "K_PLB in component BG_parameters (per_fmol)" legend_constants[49] = "K_PKACI in component BG_parameters (per_fmol)" legend_constants[50] = "K_PLB_PKACI in component BG_parameters (per_fmol)" legend_constants[51] = "K_PP1 in component BG_parameters (per_fmol)" legend_constants[52] = "K_PLBp_PP1 in component BG_parameters (per_fmol)" legend_constants[53] = "K_PLBp in component BG_parameters (per_fmol)" legend_constants[54] = "K_Ip in component BG_parameters (per_fmol)" legend_constants[55] = "K_Ip_PP1 in component BG_parameters (per_fmol)" legend_constants[56] = "K_P1_SERCA in component BG_parameters (per_fmol)" legend_constants[57] = "K_P2_SERCA in component BG_parameters (per_fmol)" legend_constants[58] = "K_P2a_SERCA in component BG_parameters (per_fmol)" legend_constants[59] = "K_P4_SERCA in component BG_parameters (per_fmol)" legend_constants[60] = "K_P5_SERCA in component BG_parameters (per_fmol)" legend_constants[61] = "K_P6_SERCA in component BG_parameters (per_fmol)" legend_constants[62] = "K_P8_SERCA in component BG_parameters (per_fmol)" legend_constants[63] = "K_P9_SERCA in component BG_parameters (per_fmol)" legend_constants[64] = "K_P10_SERCA in component BG_parameters (per_fmol)" legend_constants[65] = "K_H in component BG_parameters (per_fmol)" legend_constants[66] = "K_Cai in component BG_parameters (per_fmol)" legend_constants[67] = "K_Ca_SR in component BG_parameters (per_fmol)" legend_constants[68] = "K_MgATP in component BG_parameters (per_fmol)" legend_constants[69] = "K_MgADP in component BG_parameters (per_fmol)" legend_constants[70] = "K_P in component BG_parameters (per_fmol)" legend_constants[71] = "K_Ca_di in component BG_parameters (per_fmol)" legend_constants[72] = "K_C_RyR in component BG_parameters (per_fmol)" legend_constants[73] = "K_CI_RyR in component BG_parameters (per_fmol)" legend_constants[74] = "K_I_RyR in component BG_parameters (per_fmol)" legend_constants[75] = "K_O_RyR in component BG_parameters (per_fmol)" legend_constants[76] = "K_Cao in component BG_parameters (per_fmol)" legend_constants[77] = "K_Ki in component BG_parameters (per_fmol)" legend_constants[78] = "K_Ko in component BG_parameters (per_fmol)" legend_constants[79] = "K_000_LCC in component BG_parameters (per_fmol)" legend_constants[80] = "K_010_LCC in component BG_parameters (per_fmol)" legend_constants[81] = "K_020_LCC in component BG_parameters (per_fmol)" legend_constants[82] = "K_100_LCC in component BG_parameters (per_fmol)" legend_constants[83] = "K_110_LCC in component BG_parameters (per_fmol)" legend_constants[84] = "K_120_LCC in component BG_parameters (per_fmol)" legend_constants[85] = "K_001_LCC in component BG_parameters (per_fmol)" legend_constants[86] = "K_011_LCC in component BG_parameters (per_fmol)" legend_constants[87] = "K_021_LCC in component BG_parameters (per_fmol)" legend_constants[88] = "K_101_LCC in component BG_parameters (per_fmol)" legend_constants[89] = "K_111_LCC in component BG_parameters (per_fmol)" legend_constants[90] = "K_121_LCC in component BG_parameters (per_fmol)" legend_constants[91] = "K_TRPN in component BG_parameters (per_fmol)" legend_constants[92] = "K_Ca_TRPN in component BG_parameters (per_fmol)" legend_voi = "time in component environment (second)" legend_constants[93] = "vol_myo in component environment (pL)" legend_constants[94] = "freq in component environment (dimensionless)" legend_constants[95] = "C_m in component environment (fF)" legend_states[0] = "q_membrane in component environment (fC)" legend_algebraic[53] = "V_m in component environment (volt)" legend_algebraic[70] = "v_RyR in component RyR (fmol_per_sec)" legend_constants[96] = "F in component constants (C_per_mol)" legend_constants[97] = "q_PLB_init in component environment (fmol)" legend_constants[98] = "q_PKACI_init in component environment (fmol)" legend_constants[99] = "q_PLB_PKACI_init in component environment (fmol)" legend_constants[100] = "q_PP1_init in component environment (fmol)" legend_constants[101] = "q_PLBp_PP1_init in component environment (fmol)" legend_constants[102] = "q_PLBp_init in component environment (fmol)" legend_constants[103] = "q_Ip_init in component environment (fmol)" legend_constants[104] = "q_Ip_PP1_init in component environment (fmol)" legend_constants[105] = "q_P1_SERCA_init in component environment (fmol)" legend_constants[106] = "q_P2_SERCA_init in component environment (fmol)" legend_constants[107] = "q_P2a_SERCA_init in component environment (fmol)" legend_constants[108] = "q_P4_SERCA_init in component environment (fmol)" legend_constants[109] = "q_P5_SERCA_init in component environment (fmol)" legend_constants[110] = "q_P6_SERCA_init in component environment (fmol)" legend_constants[111] = "q_P8_SERCA_init in component environment (fmol)" legend_constants[112] = "q_P9_SERCA_init in component environment (fmol)" legend_constants[113] = "q_P10_SERCA_init in component environment (fmol)" legend_constants[114] = "q_Cai_init in component environment (fmol)" legend_constants[115] = "q_Ca_SR_init in component environment (fmol)" legend_constants[116] = "q_H_init in component environment (fmol)" legend_constants[117] = "q_P_init in component environment (fmol)" legend_constants[118] = "q_MgADP_init in component environment (fmol)" legend_constants[119] = "q_MgATP_init in component environment (fmol)" legend_constants[120] = "q_Ca_di_init in component environment (fmol)" legend_constants[121] = "q_Cao_init in component environment (fmol)" legend_constants[122] = "q_Ki_init in component environment (fmol)" legend_constants[123] = "q_Ko_init in component environment (fmol)" legend_constants[124] = "q_TRPN_init in component environment (fmol)" legend_constants[125] = "q_Ca_TRPN_init in component environment (fmol)" legend_algebraic[0] = "q_PLB in component environment (fmol)" legend_algebraic[8] = "q_PKACI in component environment (fmol)" legend_algebraic[10] = "q_PLB_PKACI in component environment (fmol)" legend_algebraic[12] = "q_PP1 in component environment (fmol)" legend_algebraic[14] = "q_PLBp_PP1 in component environment (fmol)" legend_algebraic[16] = "q_PLBp in component environment (fmol)" legend_algebraic[18] = "q_Ip in component environment (fmol)" legend_algebraic[20] = "q_Ip_PP1 in component environment (fmol)" legend_algebraic[1] = "q_P1_SERCA in component environment (fmol)" legend_algebraic[9] = "q_P2_SERCA in component environment (fmol)" legend_algebraic[11] = "q_P2a_SERCA in component environment (fmol)" legend_algebraic[13] = "q_P4_SERCA in component environment (fmol)" legend_algebraic[15] = "q_P5_SERCA in component environment (fmol)" legend_algebraic[17] = "q_P6_SERCA in component environment (fmol)" legend_algebraic[19] = "q_P8_SERCA in component environment (fmol)" legend_algebraic[21] = "q_P9_SERCA in component environment (fmol)" legend_algebraic[22] = "q_P10_SERCA in component environment (fmol)" legend_algebraic[24] = "q_H in component environment (fmol)" legend_algebraic[26] = "q_Cai in component environment (fmol)" legend_algebraic[31] = "q_MgATP in component environment (fmol)" legend_algebraic[35] = "q_MgADP in component environment (fmol)" legend_algebraic[39] = "q_P in component environment (fmol)" legend_algebraic[28] = "q_Ca_SR in component environment (fmol)" legend_algebraic[34] = "q_Ca_di in component environment (fmol)" legend_algebraic[37] = "q_Cao in component environment (fmol)" legend_algebraic[41] = "q_Ki in component environment (fmol)" legend_algebraic[45] = "q_Ko in component environment (fmol)" legend_algebraic[29] = "q_TRPN in component environment (fmol)" legend_algebraic[32] = "q_Ca_TRPN in component environment (fmol)" legend_states[1] = "q_PLB in component PLB (fmol)" legend_states[2] = "q_PKACI in component PLB (fmol)" legend_states[3] = "q_PLB_PKACI in component PLB (fmol)" legend_states[4] = "q_PP1 in component PLB (fmol)" legend_states[5] = "q_PLBp_PP1 in component PLB (fmol)" legend_states[6] = "q_PLBp in component PLB (fmol)" legend_states[7] = "q_Ip in component PLB (fmol)" legend_states[8] = "q_Ip_PP1 in component PLB (fmol)" legend_states[9] = "q_P1_SERCA in component SERCA (fmol)" legend_states[10] = "q_P2_SERCA in component SERCA (fmol)" legend_states[11] = "q_P2a_SERCA in component SERCA (fmol)" legend_states[12] = "q_P4_SERCA in component SERCA (fmol)" legend_states[13] = "q_P5_SERCA in component SERCA (fmol)" legend_states[14] = "q_P6_SERCA in component SERCA (fmol)" legend_states[15] = "q_P8_SERCA in component SERCA (fmol)" legend_states[16] = "q_P9_SERCA in component SERCA (fmol)" legend_states[17] = "q_P10_SERCA in component SERCA (fmol)" legend_states[18] = "q_H in component SERCA (fmol)" legend_states[19] = "q_Cai in component SERCA (fmol)" legend_states[20] = "q_Ca_SR in component SERCA (fmol)" legend_states[21] = "q_MgATP in component SERCA (fmol)" legend_states[22] = "q_MgADP in component SERCA (fmol)" legend_states[23] = "q_P in component SERCA (fmol)" legend_states[24] = "q_Ca_SR in component RyR (fmol)" legend_states[25] = "q_Ca_di in component RyR (fmol)" legend_states[26] = "q_Cai in component LCC (fmol)" legend_states[27] = "q_Cao in component LCC (fmol)" legend_states[28] = "q_Ki in component LCC (fmol)" legend_states[29] = "q_Ko in component LCC (fmol)" legend_states[30] = "q_Cai in component TRPN (fmol)" legend_states[31] = "q_TRPN in component TRPN (fmol)" legend_states[32] = "q_Ca_TRPN in component TRPN (fmol)" legend_algebraic[49] = "I_pulse in component environment (fA)" legend_constants[126] = "pulse_start in component environment (second)" legend_constants[127] = "pulse_end in component environment (second)" legend_constants[128] = "pulseMag in component environment (fA)" legend_constants[129] = "pulseHolding in component environment (fA)" legend_algebraic[73] = "v_E_RyR in component environment (fA)" legend_algebraic[197] = "v_E_LCC in component LCC (fA)" legend_algebraic[198] = "sum_I in component environment (fA)" legend_constants[130] = "zCa in component ion_valences (dimensionless)" legend_constants[131] = "zK in component ion_valences (dimensionless)" legend_constants[132] = "R in component constants (J_per_K_per_mol)" legend_constants[133] = "T in component constants (kelvin)" legend_constants[134] = "zNa in component ion_valences (dimensionless)" legend_constants[135] = "zCl in component ion_valences (dimensionless)" legend_algebraic[47] = "vPLBph1 in component PLB (fmol_per_sec)" legend_algebraic[51] = "vPLBph2 in component PLB (fmol_per_sec)" legend_algebraic[54] = "vPLBd1 in component PLB (fmol_per_sec)" legend_algebraic[57] = "vPLBd2 in component PLB (fmol_per_sec)" legend_algebraic[60] = "vInh in component PLB (fmol_per_sec)" legend_algebraic[23] = "mu_PLB in component PLB (J_per_mol)" legend_algebraic[25] = "mu_PKACI in component PLB (J_per_mol)" legend_algebraic[27] = "mu_PLB_PKACI in component PLB (J_per_mol)" legend_algebraic[30] = "mu_PP1 in component PLB (J_per_mol)" legend_algebraic[33] = "mu_PLBp_PP1 in component PLB (J_per_mol)" legend_algebraic[36] = "mu_PLBp in component PLB (J_per_mol)" legend_algebraic[40] = "mu_Ip in component PLB (J_per_mol)" legend_algebraic[43] = "mu_Ip_PP1 in component PLB (J_per_mol)" legend_constants[136] = "n_Cai in component SERCA (dimensionless)" legend_constants[137] = "n_Ca_SR in component SERCA (dimensionless)" legend_constants[138] = "n_H in component SERCA (dimensionless)" legend_algebraic[2] = "c_Cai in component SERCA (mM)" legend_algebraic[3] = "c_Ca_SR in component SERCA (mM)" legend_algebraic[4] = "c_H in component SERCA (mM)" legend_algebraic[5] = "c_MgADP in component SERCA (mM)" legend_algebraic[6] = "c_MgATP in component SERCA (mM)" legend_algebraic[7] = "c_P in component SERCA (mM)" legend_constants[139] = "vol_i in component SERCA (pL)" legend_constants[156] = "vol_sr in component SERCA (pL)" legend_constants[157] = "vol_isr in component SERCA (pL)" legend_algebraic[44] = "mu_Cai in component SERCA (J_per_mol)" legend_algebraic[102] = "v_Cai in component SERCA (fmol_per_sec)" legend_algebraic[48] = "mu_Ca_SR in component SERCA (J_per_mol)" legend_algebraic[137] = "v_Ca_SR in component SERCA (fmol_per_sec)" legend_algebraic[52] = "mu_H in component SERCA (J_per_mol)" legend_algebraic[152] = "v_H in component SERCA (fmol_per_sec)" legend_algebraic[55] = "mu_MgADP in component SERCA (J_per_mol)" legend_algebraic[125] = "v_MgADP in component SERCA (fmol_per_sec)" legend_algebraic[58] = "mu_MgATP in component SERCA (J_per_mol)" legend_algebraic[77] = "v_MgATP in component SERCA (fmol_per_sec)" legend_algebraic[61] = "mu_P1 in component SERCA (J_per_mol)" legend_algebraic[117] = "v_P1 in component SERCA (fmol_per_sec)" legend_algebraic[63] = "mu_P in component SERCA (J_per_mol)" legend_algebraic[118] = "v_P in component SERCA (fmol_per_sec)" legend_algebraic[69] = "mu_P2 in component SERCA (J_per_mol)" legend_algebraic[103] = "v_P2 in component SERCA (fmol_per_sec)" legend_algebraic[78] = "mu_P2a in component SERCA (J_per_mol)" legend_algebraic[95] = "v_P2a in component SERCA (fmol_per_sec)" legend_algebraic[81] = "mu_P4 in component SERCA (J_per_mol)" legend_algebraic[110] = "v_P4 in component SERCA (fmol_per_sec)" legend_algebraic[84] = "mu_P5 in component SERCA (J_per_mol)" legend_algebraic[126] = "v_P5 in component SERCA (fmol_per_sec)" legend_algebraic[119] = "mu_P6 in component SERCA (J_per_mol)" legend_algebraic[138] = "v_P6 in component SERCA (fmol_per_sec)" legend_algebraic[127] = "mu_P8 in component SERCA (J_per_mol)" legend_algebraic[146] = "v_P8 in component SERCA (fmol_per_sec)" legend_algebraic[129] = "mu_P9 in component SERCA (J_per_mol)" legend_algebraic[153] = "v_P9 in component SERCA (fmol_per_sec)" legend_algebraic[111] = "mu_P10 in component SERCA (J_per_mol)" legend_algebraic[154] = "v_P10 in component SERCA (fmol_per_sec)" legend_algebraic[67] = "Af_R1_2 in component SERCA (J_per_mol)" legend_algebraic[71] = "Ar_R1_2 in component SERCA (J_per_mol)" legend_algebraic[74] = "v_SERCA_R1_2 in component SERCA (fmol_per_sec)" legend_algebraic[87] = "Af_R5_6 in component SERCA (J_per_mol)" legend_algebraic[121] = "Ar_R5_6 in component SERCA (J_per_mol)" legend_algebraic[123] = "v_SERCA_R5_6 in component SERCA (fmol_per_sec)" legend_algebraic[89] = "Af_R2_2a in component SERCA (J_per_mol)" legend_algebraic[91] = "Ar_R2_2a in component SERCA (J_per_mol)" legend_algebraic[93] = "v_SERCA_R2_2a in component SERCA (fmol_per_sec)" legend_algebraic[96] = "Af_R2_4 in component SERCA (J_per_mol)" legend_algebraic[98] = "Ar_R2_4 in component SERCA (J_per_mol)" legend_algebraic[100] = "v_SERCA_R2_4 in component SERCA (fmol_per_sec)" legend_algebraic[104] = "Af_R4_5 in component SERCA (J_per_mol)" legend_algebraic[106] = "Ar_R4_5 in component SERCA (J_per_mol)" legend_algebraic[108] = "v_SERCA_R4_5 in component SERCA (fmol_per_sec)" legend_algebraic[131] = "Af_R6_8 in component SERCA (J_per_mol)" legend_algebraic[133] = "Ar_R6_8 in component SERCA (J_per_mol)" legend_algebraic[135] = "v_SERCA_R6_8 in component SERCA (fmol_per_sec)" legend_algebraic[139] = "Af_R8_9 in component SERCA (J_per_mol)" legend_algebraic[142] = "Ar_R8_9 in component SERCA (J_per_mol)" legend_algebraic[144] = "v_SERCA_R8_9 in component SERCA (fmol_per_sec)" legend_algebraic[147] = "Af_R9_10 in component SERCA (J_per_mol)" legend_algebraic[149] = "Ar_R9_10 in component SERCA (J_per_mol)" legend_algebraic[151] = "v_SERCA_R9_10 in component SERCA (fmol_per_sec)" legend_algebraic[113] = "Af_R10_1 in component SERCA (J_per_mol)" legend_algebraic[65] = "Ar_R10_1 in component SERCA (J_per_mol)" legend_algebraic[115] = "v_SERCA_R10_1 in component SERCA (fmol_per_sec)" legend_constants[140] = "nCa_1 in component RyR (dimensionless)" legend_constants[141] = "nCa_2 in component RyR (dimensionless)" legend_constants[142] = "q_C_RyR in component RyR (fmol)" legend_constants[143] = "q_CI_RyR in component RyR (fmol)" legend_constants[144] = "q_I_RyR in component RyR (fmol)" legend_constants[145] = "q_O_RyR in component RyR (fmol)" legend_algebraic[56] = "mu_Ca_SR in component RyR (J_per_mol)" legend_algebraic[59] = "mu_Ca_di in component RyR (J_per_mol)" legend_algebraic[85] = "v_gate_Ca_di in component RyR (fmol_per_sec)" legend_algebraic[68] = "mu_O_RyR in component RyR_gate (J_per_mol)" legend_algebraic[72] = "v_OC in component RyR_gate (fmol_per_sec)" legend_algebraic[75] = "v_CCI in component RyR_gate (fmol_per_sec)" legend_algebraic[79] = "v_CII in component RyR_gate (fmol_per_sec)" legend_algebraic[82] = "v_IO in component RyR_gate (fmol_per_sec)" legend_states[33] = "q_C_RyR in component RyR_gate (fmol)" legend_states[34] = "q_CI_RyR in component RyR_gate (fmol)" legend_states[35] = "q_I_RyR in component RyR_gate (fmol)" legend_states[36] = "q_O_RyR in component RyR_gate (fmol)" legend_algebraic[62] = "mu_C_RyR in component RyR_gate (J_per_mol)" legend_algebraic[64] = "mu_CI_RyR in component RyR_gate (J_per_mol)" legend_algebraic[66] = "mu_I_RyR in component RyR_gate (J_per_mol)" legend_constants[155] = "RT in component LCC (J_per_mol)" legend_algebraic[76] = "mu_E in component LCC (J_per_C)" legend_algebraic[90] = "mu_mem_Ca in component LCC (J_per_mol)" legend_algebraic[92] = "mu_mem_K in component LCC (J_per_mol)" legend_algebraic[165] = "v_CaGHK_i in component LCC (fmol_per_sec)" legend_algebraic[155] = "v_CaGHK_o in component LCC (fmol_per_sec)" legend_algebraic[157] = "v_KGHK_i in component LCC (fmol_per_sec)" legend_algebraic[158] = "v_KGHK_o in component LCC (fmol_per_sec)" legend_algebraic[80] = "mu_Cai in component LCC (J_per_mol)" legend_algebraic[83] = "mu_Cao in component LCC (J_per_mol)" legend_algebraic[86] = "mu_Ki in component LCC (J_per_mol)" legend_algebraic[88] = "mu_Ko in component LCC (J_per_mol)" legend_algebraic[120] = "mu_S111 in component LCC_gate (J_per_mol)" legend_algebraic[140] = "mu_S121 in component LCC_gate (J_per_mol)" legend_algebraic[122] = "Af_Ca1 in component LCC (J_per_mol)" legend_algebraic[124] = "Ar_Ca1 in component LCC (J_per_mol)" legend_algebraic[132] = "v_Ca1 in component LCC (fmol_per_sec)" legend_algebraic[150] = "v_Ca2 in component LCC (fmol_per_sec)" legend_algebraic[134] = "v_K1 in component LCC (fmol_per_sec)" legend_algebraic[156] = "v_K2 in component LCC (fmol_per_sec)" legend_algebraic[141] = "Af_Ca2 in component LCC (J_per_mol)" legend_algebraic[143] = "Ar_Ca2 in component LCC (J_per_mol)" legend_algebraic[128] = "Af_K1 in component LCC (J_per_mol)" legend_algebraic[130] = "Ar_K1 in component LCC (J_per_mol)" legend_algebraic[145] = "Af_K2 in component LCC (J_per_mol)" legend_algebraic[148] = "Ar_K2 in component LCC (J_per_mol)" legend_algebraic[159] = "v_fCa00 in component LCC_gate (fmol_per_sec)" legend_algebraic[160] = "v_fCa01 in component LCC_gate (fmol_per_sec)" legend_algebraic[161] = "v_fCa02 in component LCC_gate (fmol_per_sec)" legend_algebraic[162] = "v_fCa10 in component LCC_gate (fmol_per_sec)" legend_algebraic[163] = "v_fCa11 in component LCC_gate (fmol_per_sec)" legend_algebraic[164] = "v_fCa12 in component LCC_gate (fmol_per_sec)" legend_algebraic[194] = "v_gate in component LCC_gate (fmol_per_sec)" legend_algebraic[94] = "v_pulse in component LCC (fmol_per_sec)" legend_states[37] = "q_S000 in component LCC_gate (fmol)" legend_states[38] = "q_S010 in component LCC_gate (fmol)" legend_states[39] = "q_S020 in component LCC_gate (fmol)" legend_states[40] = "q_S100 in component LCC_gate (fmol)" legend_states[41] = "q_S110 in component LCC_gate (fmol)" legend_states[42] = "q_S120 in component LCC_gate (fmol)" legend_states[43] = "q_S001 in component LCC_gate (fmol)" legend_states[44] = "q_S011 in component LCC_gate (fmol)" legend_states[45] = "q_S021 in component LCC_gate (fmol)" legend_states[46] = "q_S101 in component LCC_gate (fmol)" legend_states[47] = "q_S111 in component LCC_gate (fmol)" legend_states[48] = "q_S121 in component LCC_gate (fmol)" legend_algebraic[97] = "mu_S000 in component LCC_gate (J_per_mol)" legend_algebraic[101] = "mu_S010 in component LCC_gate (J_per_mol)" legend_algebraic[107] = "mu_S020 in component LCC_gate (J_per_mol)" legend_algebraic[112] = "mu_S100 in component LCC_gate (J_per_mol)" legend_algebraic[116] = "mu_S110 in component LCC_gate (J_per_mol)" legend_algebraic[136] = "mu_S120 in component LCC_gate (J_per_mol)" legend_algebraic[99] = "mu_S001 in component LCC_gate (J_per_mol)" legend_algebraic[105] = "mu_S011 in component LCC_gate (J_per_mol)" legend_algebraic[109] = "mu_S021 in component LCC_gate (J_per_mol)" legend_algebraic[114] = "mu_S101 in component LCC_gate (J_per_mol)" legend_algebraic[172] = "v_f1_000 in component LCC_gate (fmol_per_sec)" legend_algebraic[173] = "v_f2_000 in component LCC_gate (fmol_per_sec)" legend_algebraic[174] = "v_f3_010 in component LCC_gate (fmol_per_sec)" legend_algebraic[176] = "v_f1_100 in component LCC_gate (fmol_per_sec)" legend_algebraic[178] = "v_f2_100 in component LCC_gate (fmol_per_sec)" legend_algebraic[180] = "v_f3_110 in component LCC_gate (fmol_per_sec)" legend_algebraic[166] = "v_d000 in component LCC_gate (fmol_per_sec)" legend_algebraic[167] = "v_d010 in component LCC_gate (fmol_per_sec)" legend_algebraic[168] = "v_d020 in component LCC_gate (fmol_per_sec)" legend_algebraic[169] = "v_d001 in component LCC_gate (fmol_per_sec)" legend_algebraic[170] = "v_d011 in component LCC_gate (fmol_per_sec)" legend_algebraic[171] = "v_d021 in component LCC_gate (fmol_per_sec)" legend_algebraic[182] = "v_f1_001 in component LCC_gate (fmol_per_sec)" legend_algebraic[185] = "v_f2_001 in component LCC_gate (fmol_per_sec)" legend_algebraic[186] = "v_f3_011 in component LCC_gate (fmol_per_sec)" legend_algebraic[188] = "v_f1_101 in component LCC_gate (fmol_per_sec)" legend_algebraic[190] = "v_f2_101 in component LCC_gate (fmol_per_sec)" legend_algebraic[192] = "v_f3_111 in component LCC_gate (fmol_per_sec)" legend_algebraic[175] = "v_S000 in component LCC_gate (fmol_per_sec)" legend_algebraic[179] = "v_S010 in component LCC_gate (fmol_per_sec)" legend_algebraic[177] = "v_S020 in component LCC_gate (fmol_per_sec)" legend_algebraic[181] = "v_S100 in component LCC_gate (fmol_per_sec)" legend_algebraic[183] = "v_S110 in component LCC_gate (fmol_per_sec)" legend_algebraic[184] = "v_S120 in component LCC_gate (fmol_per_sec)" legend_algebraic[187] = "v_S001 in component LCC_gate (fmol_per_sec)" legend_algebraic[191] = "v_S011 in component LCC_gate (fmol_per_sec)" legend_algebraic[189] = "v_S021 in component LCC_gate (fmol_per_sec)" legend_algebraic[193] = "v_S101 in component LCC_gate (fmol_per_sec)" legend_algebraic[195] = "v_S111 in component LCC_gate (fmol_per_sec)" legend_algebraic[196] = "v_S121 in component LCC_gate (fmol_per_sec)" legend_constants[146] = "zCa in component LCC_gate (dimensionless)" legend_constants[147] = "z_fd in component LCC_gate (dimensionless)" legend_constants[148] = "z_ff1 in component LCC_gate (dimensionless)" legend_constants[149] = "z_ff2 in component LCC_gate (dimensionless)" legend_constants[150] = "z_ff3 in component LCC_gate (dimensionless)" legend_constants[151] = "z_rd in component LCC_gate (dimensionless)" legend_constants[152] = "z_rf1 in component LCC_gate (dimensionless)" legend_constants[153] = "z_rf2 in component LCC_gate (dimensionless)" legend_constants[154] = "z_rf3 in component LCC_gate (dimensionless)" legend_algebraic[38] = "mu_Cai in component TRPN (J_per_mol)" legend_algebraic[42] = "mu_TRPN in component TRPN (J_per_mol)" legend_algebraic[46] = "mu_Ca_TRPN in component TRPN (J_per_mol)" legend_algebraic[50] = "v_R_TRPNCa in component TRPN (fmol_per_sec)" legend_rates[0] = "d/dt q_membrane in component environment (fC)" legend_rates[1] = "d/dt q_PLB in component PLB (fmol)" legend_rates[2] = "d/dt q_PKACI in component PLB (fmol)" legend_rates[3] = "d/dt q_PLB_PKACI in component PLB (fmol)" legend_rates[4] = "d/dt q_PP1 in component PLB (fmol)" legend_rates[5] = "d/dt q_PLBp_PP1 in component PLB (fmol)" legend_rates[6] = "d/dt q_PLBp in component PLB (fmol)" legend_rates[7] = "d/dt q_Ip in component PLB (fmol)" legend_rates[8] = "d/dt q_Ip_PP1 in component PLB (fmol)" legend_rates[19] = "d/dt q_Cai in component SERCA (fmol)" legend_rates[20] = "d/dt q_Ca_SR in component SERCA (fmol)" legend_rates[18] = "d/dt q_H in component SERCA (fmol)" legend_rates[22] = "d/dt q_MgADP in component SERCA (fmol)" legend_rates[21] = "d/dt q_MgATP in component SERCA (fmol)" legend_rates[9] = "d/dt q_P1_SERCA in component SERCA (fmol)" legend_rates[23] = "d/dt q_P in component SERCA (fmol)" legend_rates[10] = "d/dt q_P2_SERCA in component SERCA (fmol)" legend_rates[11] = "d/dt q_P2a_SERCA in component SERCA (fmol)" legend_rates[12] = "d/dt q_P4_SERCA in component SERCA (fmol)" legend_rates[13] = "d/dt q_P5_SERCA in component SERCA (fmol)" legend_rates[17] = "d/dt q_P10_SERCA in component SERCA (fmol)" legend_rates[14] = "d/dt q_P6_SERCA in component SERCA (fmol)" legend_rates[15] = "d/dt q_P8_SERCA in component SERCA (fmol)" legend_rates[16] = "d/dt q_P9_SERCA in component SERCA (fmol)" legend_rates[24] = "d/dt q_Ca_SR in component RyR (fmol)" legend_rates[25] = "d/dt q_Ca_di in component RyR (fmol)" legend_rates[36] = "d/dt q_O_RyR in component RyR_gate (fmol)" legend_rates[33] = "d/dt q_C_RyR in component RyR_gate (fmol)" legend_rates[34] = "d/dt q_CI_RyR in component RyR_gate (fmol)" legend_rates[35] = "d/dt q_I_RyR in component RyR_gate (fmol)" legend_rates[26] = "d/dt q_Cai in component LCC (fmol)" legend_rates[27] = "d/dt q_Cao in component LCC (fmol)" legend_rates[28] = "d/dt q_Ki in component LCC (fmol)" legend_rates[29] = "d/dt q_Ko in component LCC (fmol)" legend_rates[37] = "d/dt q_S000 in component LCC_gate (fmol)" legend_rates[43] = "d/dt q_S001 in component LCC_gate (fmol)" legend_rates[38] = "d/dt q_S010 in component LCC_gate (fmol)" legend_rates[44] = "d/dt q_S011 in component LCC_gate (fmol)" legend_rates[39] = "d/dt q_S020 in component LCC_gate (fmol)" legend_rates[45] = "d/dt q_S021 in component LCC_gate (fmol)" legend_rates[40] = "d/dt q_S100 in component LCC_gate (fmol)" legend_rates[46] = "d/dt q_S101 in component LCC_gate (fmol)" legend_rates[41] = "d/dt q_S110 in component LCC_gate (fmol)" legend_rates[47] = "d/dt q_S111 in component LCC_gate (fmol)" legend_rates[42] = "d/dt q_S120 in component LCC_gate (fmol)" legend_rates[48] = "d/dt q_S121 in component LCC_gate (fmol)" legend_rates[30] = "d/dt q_Cai in component TRPN (fmol)" legend_rates[31] = "d/dt q_TRPN in component TRPN (fmol)" legend_rates[32] = "d/dt q_Ca_TRPN in component TRPN (fmol)" return (legend_states, legend_algebraic, legend_voi, legend_constants) def initConsts(Gmult,igs): hackMultiplier = 1e0 constants = [0.0] * sizeConstants; states = [0.0] * sizeStates; constants[0] = 45.5263 constants[1] = 6.55904 constants[2] = 0.386674 constants[3] = 1.21269 constants[4] = 431.435 constants[5] = 1.80718e-05 constants[6] = 109628 constants[7] = 2.08432e+06 constants[8] = 109628 constants[9] = 3974.36 constants[10] = 1.35422e+06 constants[11] = 1.43841e+07 constants[12] = 1.43841e+07 constants[13] = 0.00842322 constants[14] = 1.8838e+06 constants[15] = 8.83262 constants[16] = 0.12618 constants[17] = 883.262 constants[18] = 73.6052 constants[19] = 10.2101 constants[20] = 19.0927 constants[21] = 0.0968812 constants[22] = 0.181165 constants[23] = 812.745 constants[24] = 20.0517 constants[25] = 37.4961 constants[26] = 0.190685 constants[27] = 0.0047045 constants[28] = 0.00879728 constants[29] = 1.49728 constants[30] = 7.43512 constants[31] = 0.00035129 constants[32] = 0.00174442 constants[33] = 2.74219 constants[34] = 13.617 constants[35] = 0.000643368 constants[36] = 0.0031948 constants[37] = 7703.17 constants[38] = 38252 constants[39] = 1.80731 constants[40] = 8.97463 constants[41] = 65290.3 constants[42] = 324215 constants[43] = 1610.82 constants[44] = 7998.9 constants[45] = 3012.18 constants[46] = 14957.7 constants[47] = 232.12 constants[48] = 0.00235741 constants[49] = 0.395191 constants[50] = 0.638527 constants[51] = 0.361989 constants[52] = 0.203757 constants[53] = 0.326014 constants[54] = 5.41093 constants[55] = 0.0673793 constants[56] = 853.972 constants[57] = 3243 constants[58] = 361.225 constants[59] = 6867.82 constants[60] = 0.0185784 constants[61] = 2789.64 constants[62] = 47.7378 constants[63] = 52.3431 constants[64] = 70.7485 constants[65] = 297.061 constants[66] = 0.0464875 constants[67] = 0.0992057 constants[68] = 1418.2 constants[69] = 5.12197e-06 constants[70] = 0.0175044 constants[71] = 0.0992057 constants[72] = 115.191 constants[73] = 1.15191 constants[74] = 0.00197471 constants[75] = 0.197471 constants[76] = 0.0464875 constants[77] = 0.00385122 constants[78] = 0.00385122 constants[79] = 0.0174102 constants[80] = 0.705678 constants[81] = 0.377374 constants[82] = 0.00350606 constants[83] = 0.142109 constants[84] = 0.0759954 constants[85] = 74.2064 constants[86] = 3007.77 constants[87] = 1608.46 constants[88] = 14.9436 constants[89] = 605.703 constants[90] = 323.91 constants[91] = 7.83132 constants[91] = 7.83132 constants[92] = 0.0250472 constants[93] = 34.4 constants[94] = 500 constants[95] = 1.381e5 states[0] = -13039 constants[96] = 96485 constants[97] = 4.028E+00 constants[98] = 2.234E-03 constants[99] = 1.11e-2 constants[100] = 3.382E-02 constants[101] = 1.11e-2 constants[102] = 1.11e-2 constants[103] = 1.999E-03 constants[104] = 1.11e-2 constants[105] = 0.11111 constants[106] = 0.11111 constants[107] = 0.11111 constants[108] = 0.11111 constants[109] = 0.11111 constants[110] = 0.11111 constants[111] = 0.11111 constants[112] = 0.11111 constants[113] = 0.11111 constants[114] = 6.82e-3 constants[115] = 0.641 constants[116] = 0.004028 constants[117] = 570 constants[118] = 1.3794 constants[119] = 3.8 constants[120] = 3.09e-6 constants[121] = 6.84 constants[122] = 5.51E+03 constants[123] = 2.05E+02 constants[124] = 2.57 constants[125] = 1.11e-2 states[1] = 1.11e-2 states[2] = 1.11e-2 states[3] = 1.11e-2 states[4] = 1.11e-2 states[5] = 1.11e-2 states[6] = 1.11e-2 states[7] = 1.11e-2 states[8] = 1.11e-2 states[9] = 1.11e-2 states[10] = 1.11e-2 states[11] = 1.11e-2 states[12] = 1.11e-2 states[13] = 1.11e-2 states[14] = 1.11e-2 states[15] = 1.11e-2 states[16] = 1.11e-2 states[17] = 1.11e-2 states[18] = 1.11e-2 states[19] = 1.11e-2 states[20] = 1.11e-2 states[21] = 1.11e-2 states[22] = 1.11e-2 states[23] = 1.11e-2 states[24] = 1.11e-2 states[25] = 1.11e-2 states[26] = 1.11e-2 states[27] = 1.11e-2 states[28] = 1.11e-2 states[29] = 1.11e-2 states[30] = 1.11e-2 states[31] = 1.11e-2 states[32] = 1.11e-2 constants[126] = 1e-3 constants[127] = 2e-3 constants[128] = 4e7 constants[129] = 0 constants[130] = 2 constants[131] = 1 constants[132] = 8.31 constants[133] = 310 constants[134] = 1 constants[135] = -1 constants[136] = 2 constants[137] = 2 constants[138] = 2 constants[139] = 34.0 constants[140] = 1 constants[141] = 2 constants[142] = 1.11e-2 constants[143] = 1.11e-2 constants[144] = 1.11e-2 constants[145] = 1.11e-2 states[33] = 1.11e-2 states[34] = 1.11e-2 states[35] = 1.11e-2 states[36] = 1.11e-2 states[37] = 3.6989e-07 states[38] = 7.3239e-05 states[39] = 3.6989e-07 states[40] = 3.7363e-09 states[41] = 7.3979e-07 states[42] = 3.7363e-09 states[43] = 4.1099e-08 states[44] = 8.1377e-06 states[45] = 4.1099e-08 states[46] = 4.1514e-10 states[47] = 8.2199e-08 states[48] = 4.1514e-10 constants[146] = 2 constants[147] = 2.1404 constants[148] = -1.1495 constants[149] = 0.72162 constants[150] = 4.2933 constants[151] = -2.1404 constants[152] = 1.8993 constants[153] = -0.52288 constants[154] = 0 constants[155] = constants[132]*constants[133] constants[156] = constants[139]*0.0350000 constants[158] = 0.000000 constants[159] = 0.000000 constants[160] = 0.000000 constants[161] = 0.000000 constants[162] = 0.000000 constants[163] = 0.000000 constants[157] = constants[139]+constants[156] # hack increase constants specified by igs by a multiplier Gmult constants = [c*Gmult if i in igs else c for i, c in enumerate(constants)] states = [s*hackMultiplier for s in states] return (states, constants) def computeRates(voi, states, constants): print(voi) rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic rates[19] = constants[158] rates[20] = constants[159] rates[18] = constants[160] rates[22] = constants[161] rates[21] = constants[162] rates[23] = constants[163] algebraic[26] = constants[114]+states[19]+states[26]+states[30] algebraic[38] = constants[132]*constants[133]*log(constants[66]*algebraic[26]) algebraic[29] = constants[124]+states[31] algebraic[42] = constants[132]*constants[133]*log(constants[91]*algebraic[29]) algebraic[32] = constants[125]+states[32] algebraic[46] = constants[132]*constants[133]*log(constants[92]*algebraic[32]) algebraic[50] = constants[47]*(exp(algebraic[38]+algebraic[42])/(constants[132]*constants[133])-exp(algebraic[46])/(constants[132]*constants[133])) rates[30] = -algebraic[50] rates[31] = -algebraic[50] rates[32] = algebraic[50] algebraic[0] = constants[97]+states[1] algebraic[23] = constants[132]*constants[133]*log(constants[48]*algebraic[0]) algebraic[8] = constants[98]+states[2] algebraic[25] = constants[132]*constants[133]*log(constants[49]*algebraic[8]) algebraic[10] = constants[99]+states[3] algebraic[27] = constants[132]*constants[133]*log(constants[50]*algebraic[10]) algebraic[47] = constants[0]*(exp((algebraic[23]+algebraic[25])/(constants[132]*constants[133]))-exp(algebraic[27]/(constants[132]*constants[133]))) algebraic[16] = constants[102]+states[6] algebraic[36] = constants[132]*constants[133]*log(constants[53]*algebraic[16]) algebraic[51] = constants[1]*(exp(algebraic[27]/(constants[132]*constants[133]))-exp((algebraic[36]+algebraic[25])/(constants[132]*constants[133]))) rates[2] = algebraic[51]-algebraic[47] rates[3] = algebraic[47]-algebraic[51] algebraic[12] = constants[100]+states[4] algebraic[30] = constants[132]*constants[133]*log(constants[51]*algebraic[12]) algebraic[14] = constants[101]+states[5] algebraic[33] = constants[132]*constants[133]*log(constants[52]*algebraic[14]) algebraic[54] = constants[2]*(exp((algebraic[36]+algebraic[30])/(constants[132]*constants[133]))-exp(algebraic[33]/(constants[132]*constants[133]))) rates[6] = algebraic[51]-algebraic[54] algebraic[57] = constants[3]*(exp(algebraic[33]/(constants[132]*constants[133]))-exp((algebraic[23]+algebraic[30])/(constants[132]*constants[133]))) rates[1] = algebraic[57]-algebraic[47] rates[5] = algebraic[54]-algebraic[57] algebraic[18] = constants[103]+states[7] algebraic[40] = constants[132]*constants[133]*log(constants[54]*algebraic[18]) algebraic[20] = constants[104]+states[8] algebraic[43] = constants[132]*constants[133]*log(constants[55]*algebraic[20]) algebraic[60] = constants[4]*(exp((algebraic[30]+algebraic[40])/(constants[132]*constants[133]))-exp(algebraic[43]/(constants[132]*constants[133]))) rates[4] = (algebraic[57]-algebraic[54])-algebraic[60] rates[7] = -algebraic[60] rates[8] = algebraic[60] algebraic[53] = states[0]/constants[95] algebraic[28] = constants[115]+states[20]+states[24] algebraic[56] = constants[132]*constants[133]*log(constants[67]*algebraic[28]) algebraic[34] = constants[120]+states[25] algebraic[59] = constants[132]*constants[133]*log(constants[71]*algebraic[34]) algebraic[68] = constants[132]*constants[133]*log(constants[75]*states[36]) algebraic[70] = constants[14]*exp(algebraic[68]/(constants[132]*constants[133]))*(exp((algebraic[56]+constants[130]*constants[96]*algebraic[53])/(constants[132]*constants[133]))-exp(algebraic[59]/(constants[132]*constants[133]))) rates[24] = -algebraic[70] algebraic[62] = constants[132]*constants[133]*log(constants[72]*states[33]) algebraic[72] = constants[15]*(exp((algebraic[68]+constants[141]*algebraic[59])/(constants[132]*constants[133]))-exp(algebraic[62]/(constants[132]*constants[133]))) algebraic[64] = constants[132]*constants[133]*log(constants[73]*states[34]) algebraic[75] = constants[16]*(exp((algebraic[62]+constants[140]*algebraic[59])/(constants[132]*constants[133]))-exp(algebraic[64]/(constants[132]*constants[133]))) rates[33] = algebraic[72]-algebraic[75] algebraic[66] = constants[132]*constants[133]*log(constants[74]*states[35]) algebraic[79] = constants[17]*(exp((algebraic[64]+constants[141]*algebraic[59])/(constants[132]*constants[133]))-exp(algebraic[66]/(constants[132]*constants[133]))) rates[34] = algebraic[75]-algebraic[79] algebraic[82] = constants[18]*(exp((algebraic[66]+constants[140]*algebraic[59])/(constants[132]*constants[133]))-exp(algebraic[68]/(constants[132]*constants[133]))) rates[36] = algebraic[82]-algebraic[72] rates[35] = algebraic[79]-algebraic[82] algebraic[85] = ((-constants[141]*algebraic[72]-constants[140]*algebraic[75])-constants[141]*algebraic[79])-constants[140]*algebraic[82] rates[25] = algebraic[70]+algebraic[85] algebraic[24] = constants[116]+states[18] algebraic[52] = constants[132]*constants[133]*log(constants[65]*algebraic[24]) algebraic[9] = constants[106]+states[10] algebraic[69] = constants[132]*constants[133]*log(constants[57]*algebraic[9]) algebraic[89] = algebraic[69]+algebraic[52] algebraic[11] = constants[107]+states[11] algebraic[78] = constants[132]*constants[133]*log(constants[58]*algebraic[11]) algebraic[91] = algebraic[78] algebraic[93] = constants[7]*(exp(algebraic[89]/(constants[132]*constants[133]))-exp(algebraic[91]/(constants[132]*constants[133]))) algebraic[95] = algebraic[93] rates[11] = algebraic[95] algebraic[31] = constants[119]+states[21] algebraic[58] = constants[132]*constants[133]*log(constants[68]*algebraic[31]) algebraic[1] = constants[105]+states[9] algebraic[61] = constants[132]*constants[133]*log(constants[56]*algebraic[1]) algebraic[67] = algebraic[61]+algebraic[58] algebraic[71] = algebraic[69] algebraic[74] = constants[5]*(exp(algebraic[67]/(constants[132]*constants[133]))-exp(algebraic[71]/(constants[132]*constants[133]))) algebraic[44] = constants[132]*constants[133]*log(constants[66]*algebraic[26]) algebraic[96] = algebraic[69]+constants[136]*algebraic[44] algebraic[13] = constants[108]+states[12] algebraic[81] = constants[132]*constants[133]*log(constants[59]*algebraic[13]) algebraic[98] = algebraic[81] algebraic[100] = constants[6]*(exp(algebraic[96]/(constants[132]*constants[133]))-exp(algebraic[98]/(constants[132]*constants[133]))) algebraic[103] = algebraic[74]-algebraic[100] rates[10] = algebraic[103] algebraic[104] = algebraic[81] algebraic[15] = constants[109]+states[13] algebraic[84] = constants[132]*constants[133]*log(constants[60]*algebraic[15]) algebraic[106] = algebraic[84]+constants[138]*algebraic[52] algebraic[108] = constants[8]*(exp(algebraic[104]/(constants[132]*constants[133]))-exp(algebraic[106]/(constants[132]*constants[133]))) algebraic[110] = algebraic[100]-algebraic[108] rates[12] = algebraic[110] algebraic[22] = constants[113]+states[17] algebraic[111] = constants[132]*constants[133]*log(constants[64]*algebraic[22]) algebraic[113] = algebraic[111] algebraic[39] = constants[117]+states[23] algebraic[63] = constants[132]*constants[133]*log(constants[70]*algebraic[39]) algebraic[65] = algebraic[61]+algebraic[63] algebraic[115] = constants[13]*(exp(algebraic[113]/(constants[132]*constants[133]))-exp(algebraic[65]/(constants[132]*constants[133]))) algebraic[117] = algebraic[115]-algebraic[74] rates[9] = algebraic[117] algebraic[87] = algebraic[84] algebraic[35] = constants[118]+states[22] algebraic[55] = constants[132]*constants[133]*log(constants[69]*algebraic[35]) algebraic[17] = constants[110]+states[14] algebraic[119] = constants[132]*constants[133]*log(constants[61]*algebraic[17]) algebraic[121] = algebraic[55]+algebraic[119] algebraic[123] = constants[9]*(exp(algebraic[87]/(constants[132]*constants[133]))-exp(algebraic[121]/(constants[132]*constants[133]))) algebraic[126] = algebraic[108]-algebraic[123] rates[13] = algebraic[126] algebraic[131] = algebraic[119] algebraic[48] = constants[132]*constants[133]*log(constants[67]*algebraic[28]) algebraic[19] = constants[111]+states[15] algebraic[127] = constants[132]*constants[133]*log(constants[62]*algebraic[19]) algebraic[133] = algebraic[127]+constants[137]*algebraic[48] algebraic[135] = constants[10]*(exp(algebraic[131]/(constants[132]*constants[133]))-exp(algebraic[133]/(constants[132]*constants[133]))) algebraic[138] = algebraic[123]-algebraic[135] rates[14] = algebraic[138] algebraic[139] = algebraic[127]+constants[138]*algebraic[52] algebraic[21] = constants[112]+states[16] algebraic[129] = constants[132]*constants[133]*log(constants[63]*algebraic[21]) algebraic[142] = algebraic[129] algebraic[144] = constants[11]*(exp(algebraic[139]/(constants[132]*constants[133]))-exp(algebraic[142]/(constants[132]*constants[133]))) algebraic[146] = algebraic[135]-algebraic[144] rates[15] = algebraic[146] algebraic[147] = algebraic[129] algebraic[149] = algebraic[52]+algebraic[111] algebraic[151] = constants[12]*(exp(algebraic[147]/(constants[132]*constants[133]))-exp(algebraic[149]/(constants[132]*constants[133]))) algebraic[154] = algebraic[151]-algebraic[115] rates[17] = algebraic[154] algebraic[153] = algebraic[144]-algebraic[151] rates[16] = algebraic[153] algebraic[76] = states[0]/constants[95] algebraic[90] = constants[130]*constants[96]*algebraic[76] algebraic[80] = constants[155]*log(constants[66]*algebraic[26]) algebraic[120] = constants[155]*log(constants[89]*states[47]) algebraic[122] = algebraic[80]+algebraic[90]+algebraic[120] algebraic[37] = constants[121]+states[27] algebraic[83] = constants[155]*log(constants[76]*algebraic[37]) algebraic[124] = algebraic[83]+algebraic[120] algebraic[132] = custom_piecewise([equal(algebraic[90] , 0.000000), constants[19]*(exp(algebraic[122]/constants[155])-exp(algebraic[124]/constants[155])) , True, (((constants[19]*algebraic[90])/constants[155])/(exp(algebraic[90]/constants[155])-1.00000))*(exp(algebraic[122]/constants[155])-exp(algebraic[124]/constants[155]))]) algebraic[140] = constants[155]*log(constants[90]*states[48]) algebraic[141] = algebraic[80]+algebraic[90]+algebraic[140] algebraic[143] = algebraic[83]+algebraic[140] algebraic[150] = custom_piecewise([equal(algebraic[90] , 0.000000), constants[20]*(exp(algebraic[141]/constants[155])-exp(algebraic[143]/constants[155])) , True, (((constants[20]*algebraic[90])/constants[155])/(exp(algebraic[90]/constants[155])-1.00000))*(exp(algebraic[141]/constants[155])-exp(algebraic[143]/constants[155]))]) algebraic[155] = algebraic[150]+algebraic[132] algebraic[49] = custom_piecewise([greater(voi , constants[126]) & less(voi , constants[127]), constants[128] , True, constants[129]]) algebraic[94] = algebraic[49]/(constants[130]*constants[96]) rates[27] = algebraic[155]-algebraic[94] algebraic[92] = constants[131]*constants[96]*algebraic[76] algebraic[41] = constants[122]+states[28] algebraic[86] = constants[155]*log(constants[77]*algebraic[41]) algebraic[128] = algebraic[86]+algebraic[92]+algebraic[120] algebraic[45] = constants[123]+states[29] algebraic[88] = constants[155]*log(constants[78]*algebraic[45]) algebraic[130] = algebraic[88]+algebraic[120] algebraic[134] = custom_piecewise([equal(algebraic[92] , 0.000000), constants[21]*(exp(algebraic[128]/constants[155])-exp(algebraic[130]/constants[155])) , True, (((constants[21]*algebraic[92])/constants[155])/(exp(algebraic[92]/constants[155])-1.00000))*(exp(algebraic[128]/constants[155])-exp(algebraic[130]/constants[155]))]) algebraic[145] = algebraic[86]+algebraic[92]+algebraic[140] algebraic[148] = algebraic[88]+algebraic[140] algebraic[156] = custom_piecewise([equal(algebraic[92] , 0.000000), constants[22]*(exp(algebraic[145]/constants[155])-exp(algebraic[148]/constants[155])) , True, (((constants[22]*algebraic[92])/constants[155])/(exp(algebraic[92]/constants[155])-1.00000))*(exp(algebraic[145]/constants[155])-exp(algebraic[148]/constants[155]))]) algebraic[157] = -algebraic[134]-algebraic[156] rates[28] = algebraic[157] algebraic[158] = algebraic[156]+algebraic[134] rates[29] = algebraic[158] algebraic[97] = constants[155]*log(constants[79]*states[37]) algebraic[99] = constants[155]*log(constants[85]*states[43]) algebraic[159] = constants[41]*(exp(algebraic[97]/constants[155])-exp((algebraic[99]+constants[146]*algebraic[80])/constants[155])) algebraic[101] = constants[155]*log(constants[80]*states[38]) algebraic[105] = constants[155]*log(constants[86]*states[44]) algebraic[160] = constants[43]*(exp(algebraic[101]/constants[155])-exp((algebraic[105]+constants[146]*algebraic[80])/constants[155])) algebraic[107] = constants[155]*log(constants[81]*states[39]) algebraic[109] = constants[155]*log(constants[87]*states[45]) algebraic[161] = constants[45]*(exp(algebraic[107]/constants[155])-exp((algebraic[109]+constants[146]*algebraic[80])/constants[155])) algebraic[112] = constants[155]*log(constants[82]*states[40]) algebraic[114] = constants[155]*log(constants[88]*states[46]) algebraic[162] = constants[42]*(exp(algebraic[112]/constants[155])-exp((algebraic[114]+constants[146]*algebraic[80])/constants[155])) algebraic[116] = constants[155]*log(constants[83]*states[41]) algebraic[163] = constants[44]*(exp(algebraic[116]/constants[155])-exp((algebraic[120]+constants[146]*algebraic[80])/constants[155])) algebraic[136] = constants[155]*log(constants[84]*states[42]) algebraic[164] = constants[46]*(exp(algebraic[136]/constants[155])-exp((algebraic[140]+constants[146]*algebraic[80])/constants[155])) algebraic[165] = (constants[130]*(algebraic[162]+algebraic[163]+algebraic[161]+algebraic[160]+algebraic[159]+algebraic[164])-algebraic[132])-algebraic[150] rates[26] = algebraic[165]+algebraic[94] algebraic[172] = constants[29]*(exp((algebraic[97]+constants[148]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[101]+constants[152]*constants[96]*algebraic[76])/constants[155])) algebraic[173] = constants[33]*(exp((algebraic[97]+constants[149]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[107]+constants[153]*constants[96]*algebraic[76])/constants[155])) algebraic[166] = constants[23]*(exp((algebraic[97]+constants[147]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[112]+constants[151]*constants[96]*algebraic[76])/constants[155])) algebraic[175] = ((-algebraic[166]-algebraic[172])-algebraic[173])-algebraic[159] rates[37] = algebraic[175] algebraic[174] = constants[37]*(exp((algebraic[101]+constants[150]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[107]+constants[154]*constants[96]*algebraic[76])/constants[155])) algebraic[168] = constants[25]*(exp((algebraic[107]+constants[147]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[136]+constants[151]*constants[96]*algebraic[76])/constants[155])) algebraic[177] = (-algebraic[168]+algebraic[174]+algebraic[173])-algebraic[161] rates[39] = algebraic[177] algebraic[176] = constants[30]*(exp((algebraic[112]+constants[148]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[116]+constants[152]*constants[96]*algebraic[76])/constants[155])) algebraic[167] = constants[24]*(exp((algebraic[101]+constants[147]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[116]+constants[151]*constants[96]*algebraic[76])/constants[155])) algebraic[179] = ((-algebraic[167]-algebraic[174])+algebraic[176])-algebraic[160] rates[38] = algebraic[179] algebraic[178] = constants[34]*(exp((algebraic[112]+constants[149]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[136]+constants[153]*constants[96]*algebraic[76])/constants[155])) algebraic[181] = ((-algebraic[178]+algebraic[166])-algebraic[176])-algebraic[162] rates[40] = algebraic[181] algebraic[180] = constants[38]*(exp((algebraic[116]+constants[150]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[136]+constants[154]*constants[96]*algebraic[76])/constants[155])) algebraic[183] = ((algebraic[176]+algebraic[167])-algebraic[163])-algebraic[180] rates[41] = algebraic[183] algebraic[184] = ((algebraic[178]+algebraic[168])-algebraic[164])+algebraic[180] rates[42] = algebraic[184] algebraic[169] = constants[26]*(exp((algebraic[99]+constants[147]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[114]+constants[151]*constants[96]*algebraic[76])/constants[155])) algebraic[182] = constants[31]*(exp((algebraic[99]+constants[148]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[105]+constants[152]*constants[96]*algebraic[76])/constants[155])) algebraic[185] = constants[35]*(exp((algebraic[99]+constants[149]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[109]+constants[153]*constants[96]*algebraic[76])/constants[155])) algebraic[187] = ((-algebraic[169]-algebraic[182])-algebraic[185])+algebraic[159] rates[43] = algebraic[187] algebraic[171] = constants[28]*(exp((algebraic[109]+constants[147]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[140]+constants[151]*constants[96]*algebraic[76])/constants[155])) algebraic[186] = constants[39]*(exp((algebraic[105]+constants[150]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[109]+constants[154]*constants[96]*algebraic[76])/constants[155])) algebraic[189] = -algebraic[171]+algebraic[186]+algebraic[185]+algebraic[161] rates[45] = algebraic[189] algebraic[170] = constants[27]*(exp((algebraic[105]+constants[147]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[120]+constants[151]*constants[96]*algebraic[76])/constants[155])) algebraic[188] = constants[32]*(exp((algebraic[114]+constants[148]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[120]+constants[152]*constants[96]*algebraic[76])/constants[155])) algebraic[191] = (-algebraic[170]-algebraic[186])+algebraic[188]+algebraic[160] rates[44] = algebraic[191] algebraic[190] = constants[36]*(exp((algebraic[114]+constants[149]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[140]+constants[153]*constants[96]*algebraic[76])/constants[155])) algebraic[193] = ((-algebraic[190]+algebraic[169])-algebraic[188])+algebraic[162] rates[46] = algebraic[193] algebraic[192] = constants[40]*(exp((algebraic[120]+constants[150]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[140]+constants[154]*constants[96]*algebraic[76])/constants[155])) algebraic[195] = (algebraic[188]+algebraic[170]+algebraic[163])-algebraic[192] rates[47] = algebraic[195] algebraic[196] = algebraic[190]+algebraic[171]+algebraic[164]+algebraic[192] rates[48] = algebraic[196] algebraic[73] = constants[130]*algebraic[70]*constants[96] algebraic[194] = algebraic[185]*(constants[153]-constants[149])+algebraic[186]*(constants[154]-constants[150])+algebraic[188]*(constants[152]-constants[148])+algebraic[190]*(constants[153]-constants[149])+algebraic[192]*(constants[154]-constants[150])+algebraic[182]*(constants[152]-constants[148])+algebraic[171]*(constants[151]-constants[147])+algebraic[170]*(constants[151]-constants[147])+algebraic[169]*(constants[151]-constants[147])+algebraic[180]*(constants[154]-constants[150])+algebraic[178]*(constants[153]-constants[149])+algebraic[176]*(constants[152]-constants[148])+algebraic[174]*(constants[154]-constants[150])+algebraic[173]*(constants[153]-constants[149])+algebraic[172]*(constants[152]-constants[148])+algebraic[168]*(constants[151]-constants[147])+algebraic[167]*(constants[151]-constants[147])+algebraic[166]*(constants[151]-constants[147]) algebraic[197] = constants[96]*((algebraic[194]-constants[130]*(algebraic[132]+algebraic[150]))-constants[131]*(algebraic[134]+algebraic[156])) algebraic[198] = (-algebraic[73]-algebraic[197])+algebraic[49] rates[0] = algebraic[198] return(rates) def computeAlgebraic(constants, states, voi): algebraic = array([[0.0] * len(voi)] * sizeAlgebraic) states = array(states) voi = array(voi) algebraic[26] = constants[114]+states[19]+states[26]+states[30] algebraic[38] = constants[132]*constants[133]*log(constants[66]*algebraic[26]) algebraic[29] = constants[124]+states[31] algebraic[42] = constants[132]*constants[133]*log(constants[91]*algebraic[29]) algebraic[32] = constants[125]+states[32] algebraic[46] = constants[132]*constants[133]*log(constants[92]*algebraic[32]) algebraic[50] = constants[47]*(exp(algebraic[38]+algebraic[42])/(constants[132]*constants[133])-exp(algebraic[46])/(constants[132]*constants[133])) algebraic[0] = constants[97]+states[1] algebraic[23] = constants[132]*constants[133]*log(constants[48]*algebraic[0]) algebraic[8] = constants[98]+states[2] algebraic[25] = constants[132]*constants[133]*log(constants[49]*algebraic[8]) algebraic[10] = constants[99]+states[3] algebraic[27] = constants[132]*constants[133]*log(constants[50]*algebraic[10]) algebraic[47] = constants[0]*(exp((algebraic[23]+algebraic[25])/(constants[132]*constants[133]))-exp(algebraic[27]/(constants[132]*constants[133]))) algebraic[16] = constants[102]+states[6] algebraic[36] = constants[132]*constants[133]*log(constants[53]*algebraic[16]) algebraic[51] = constants[1]*(exp(algebraic[27]/(constants[132]*constants[133]))-exp((algebraic[36]+algebraic[25])/(constants[132]*constants[133]))) algebraic[12] = constants[100]+states[4] algebraic[30] = constants[132]*constants[133]*log(constants[51]*algebraic[12]) algebraic[14] = constants[101]+states[5] algebraic[33] = constants[132]*constants[133]*log(constants[52]*algebraic[14]) algebraic[54] = constants[2]*(exp((algebraic[36]+algebraic[30])/(constants[132]*constants[133]))-exp(algebraic[33]/(constants[132]*constants[133]))) algebraic[57] = constants[3]*(exp(algebraic[33]/(constants[132]*constants[133]))-exp((algebraic[23]+algebraic[30])/(constants[132]*constants[133]))) algebraic[18] = constants[103]+states[7] algebraic[40] = constants[132]*constants[133]*log(constants[54]*algebraic[18]) algebraic[20] = constants[104]+states[8] algebraic[43] = constants[132]*constants[133]*log(constants[55]*algebraic[20]) algebraic[60] = constants[4]*(exp((algebraic[30]+algebraic[40])/(constants[132]*constants[133]))-exp(algebraic[43]/(constants[132]*constants[133]))) algebraic[53] = states[0]/constants[95] algebraic[28] = constants[115]+states[20]+states[24] algebraic[56] = constants[132]*constants[133]*log(constants[67]*algebraic[28]) algebraic[34] = constants[120]+states[25] algebraic[59] = constants[132]*constants[133]*log(constants[71]*algebraic[34]) if len(states[0]) == 1: if states[36]<=0: j = 10 algebraic[68] = constants[132]*constants[133]*log(constants[75]*states[36]) algebraic[70] = constants[14]*exp(algebraic[68]/(constants[132]*constants[133]))*(exp((algebraic[56]+constants[130]*constants[96]*algebraic[53])/(constants[132]*constants[133]))-exp(algebraic[59]/(constants[132]*constants[133]))) algebraic[62] = constants[132]*constants[133]*log(constants[72]*states[33]) algebraic[72] = constants[15]*(exp((algebraic[68]+constants[141]*algebraic[59])/(constants[132]*constants[133]))-exp(algebraic[62]/(constants[132]*constants[133]))) algebraic[64] = constants[132]*constants[133]*log(constants[73]*states[34]) algebraic[75] = constants[16]*(exp((algebraic[62]+constants[140]*algebraic[59])/(constants[132]*constants[133]))-exp(algebraic[64]/(constants[132]*constants[133]))) algebraic[66] = constants[132]*constants[133]*log(constants[74]*states[35]) algebraic[79] = constants[17]*(exp((algebraic[64]+constants[141]*algebraic[59])/(constants[132]*constants[133]))-exp(algebraic[66]/(constants[132]*constants[133]))) algebraic[82] = constants[18]*(exp((algebraic[66]+constants[140]*algebraic[59])/(constants[132]*constants[133]))-exp(algebraic[68]/(constants[132]*constants[133]))) algebraic[85] = ((-constants[141]*algebraic[72]-constants[140]*algebraic[75])-constants[141]*algebraic[79])-constants[140]*algebraic[82] algebraic[24] = constants[116]+states[18] algebraic[52] = constants[132]*constants[133]*log(constants[65]*algebraic[24]) algebraic[9] = constants[106]+states[10] algebraic[69] = constants[132]*constants[133]*log(constants[57]*algebraic[9]) algebraic[89] = algebraic[69]+algebraic[52] algebraic[11] = constants[107]+states[11] algebraic[78] = constants[132]*constants[133]*log(constants[58]*algebraic[11]) algebraic[91] = algebraic[78] algebraic[93] = constants[7]*(exp(algebraic[89]/(constants[132]*constants[133]))-exp(algebraic[91]/(constants[132]*constants[133]))) algebraic[95] = algebraic[93] algebraic[31] = constants[119]+states[21] algebraic[58] = constants[132]*constants[133]*log(constants[68]*algebraic[31]) algebraic[1] = constants[105]+states[9] algebraic[61] = constants[132]*constants[133]*log(constants[56]*algebraic[1]) algebraic[67] = algebraic[61]+algebraic[58] algebraic[71] = algebraic[69] algebraic[74] = constants[5]*(exp(algebraic[67]/(constants[132]*constants[133]))-exp(algebraic[71]/(constants[132]*constants[133]))) algebraic[44] = constants[132]*constants[133]*log(constants[66]*algebraic[26]) algebraic[96] = algebraic[69]+constants[136]*algebraic[44] algebraic[13] = constants[108]+states[12] algebraic[81] = constants[132]*constants[133]*log(constants[59]*algebraic[13]) algebraic[98] = algebraic[81] algebraic[100] = constants[6]*(exp(algebraic[96]/(constants[132]*constants[133]))-exp(algebraic[98]/(constants[132]*constants[133]))) algebraic[103] = algebraic[74]-algebraic[100] algebraic[104] = algebraic[81] algebraic[15] = constants[109]+states[13] algebraic[84] = constants[132]*constants[133]*log(constants[60]*algebraic[15]) algebraic[106] = algebraic[84]+constants[138]*algebraic[52] algebraic[108] = constants[8]*(exp(algebraic[104]/(constants[132]*constants[133]))-exp(algebraic[106]/(constants[132]*constants[133]))) algebraic[110] = algebraic[100]-algebraic[108] algebraic[22] = constants[113]+states[17] algebraic[111] = constants[132]*constants[133]*log(constants[64]*algebraic[22]) algebraic[113] = algebraic[111] algebraic[39] = constants[117]+states[23] algebraic[63] = constants[132]*constants[133]*log(constants[70]*algebraic[39]) algebraic[65] = algebraic[61]+algebraic[63] algebraic[115] = constants[13]*(exp(algebraic[113]/(constants[132]*constants[133]))-exp(algebraic[65]/(constants[132]*constants[133]))) algebraic[117] = algebraic[115]-algebraic[74] algebraic[87] = algebraic[84] algebraic[35] = constants[118]+states[22] algebraic[55] = constants[132]*constants[133]*log(constants[69]*algebraic[35]) algebraic[17] = constants[110]+states[14] algebraic[119] = constants[132]*constants[133]*log(constants[61]*algebraic[17]) algebraic[121] = algebraic[55]+algebraic[119] algebraic[123] = constants[9]*(exp(algebraic[87]/(constants[132]*constants[133]))-exp(algebraic[121]/(constants[132]*constants[133]))) algebraic[126] = algebraic[108]-algebraic[123] algebraic[131] = algebraic[119] algebraic[48] = constants[132]*constants[133]*log(constants[67]*algebraic[28]) algebraic[19] = constants[111]+states[15] algebraic[127] = constants[132]*constants[133]*log(constants[62]*algebraic[19]) algebraic[133] = algebraic[127]+constants[137]*algebraic[48] algebraic[135] = constants[10]*(exp(algebraic[131]/(constants[132]*constants[133]))-exp(algebraic[133]/(constants[132]*constants[133]))) algebraic[138] = algebraic[123]-algebraic[135] algebraic[139] = algebraic[127]+constants[138]*algebraic[52] algebraic[21] = constants[112]+states[16] algebraic[129] = constants[132]*constants[133]*log(constants[63]*algebraic[21]) algebraic[142] = algebraic[129] algebraic[144] = constants[11]*(exp(algebraic[139]/(constants[132]*constants[133]))-exp(algebraic[142]/(constants[132]*constants[133]))) algebraic[146] = algebraic[135]-algebraic[144] algebraic[147] = algebraic[129] algebraic[149] = algebraic[52]+algebraic[111] algebraic[151] = constants[12]*(exp(algebraic[147]/(constants[132]*constants[133]))-exp(algebraic[149]/(constants[132]*constants[133]))) algebraic[154] = algebraic[151]-algebraic[115] algebraic[153] = algebraic[144]-algebraic[151] algebraic[76] = states[0]/constants[95] algebraic[90] = constants[130]*constants[96]*algebraic[76] algebraic[80] = constants[155]*log(constants[66]*algebraic[26]) algebraic[120] = constants[155]*log(constants[89]*states[47]) algebraic[122] = algebraic[80]+algebraic[90]+algebraic[120] algebraic[37] = constants[121]+states[27] algebraic[83] = constants[155]*log(constants[76]*algebraic[37]) algebraic[124] = algebraic[83]+algebraic[120] algebraic[132] = custom_piecewise([equal(algebraic[90] , 0.000000), constants[19]*(exp(algebraic[122]/constants[155])-exp(algebraic[124]/constants[155])) , True, (((constants[19]*algebraic[90])/constants[155])/(exp(algebraic[90]/constants[155])-1.00000))*(exp(algebraic[122]/constants[155])-exp(algebraic[124]/constants[155]))]) algebraic[140] = constants[155]*log(constants[90]*states[48]) algebraic[141] = algebraic[80]+algebraic[90]+algebraic[140] algebraic[143] = algebraic[83]+algebraic[140] algebraic[150] = custom_piecewise([equal(algebraic[90] , 0.000000), constants[20]*(exp(algebraic[141]/constants[155])-exp(algebraic[143]/constants[155])) , True, (((constants[20]*algebraic[90])/constants[155])/(exp(algebraic[90]/constants[155])-1.00000))*(exp(algebraic[141]/constants[155])-exp(algebraic[143]/constants[155]))]) algebraic[155] = algebraic[150]+algebraic[132] algebraic[49] = custom_piecewise([greater(voi , constants[126]) & less(voi , constants[127]), constants[128] , True, constants[129]]) algebraic[94] = algebraic[49]/(constants[130]*constants[96]) algebraic[92] = constants[131]*constants[96]*algebraic[76] algebraic[41] = constants[122]+states[28] algebraic[86] = constants[155]*log(constants[77]*algebraic[41]) algebraic[128] = algebraic[86]+algebraic[92]+algebraic[120] algebraic[45] = constants[123]+states[29] algebraic[88] = constants[155]*log(constants[78]*algebraic[45]) algebraic[130] = algebraic[88]+algebraic[120] algebraic[134] = custom_piecewise([equal(algebraic[92] , 0.000000), constants[21]*(exp(algebraic[128]/constants[155])-exp(algebraic[130]/constants[155])) , True, (((constants[21]*algebraic[92])/constants[155])/(exp(algebraic[92]/constants[155])-1.00000))*(exp(algebraic[128]/constants[155])-exp(algebraic[130]/constants[155]))]) algebraic[145] = algebraic[86]+algebraic[92]+algebraic[140] algebraic[148] = algebraic[88]+algebraic[140] algebraic[156] = custom_piecewise([equal(algebraic[92] , 0.000000), constants[22]*(exp(algebraic[145]/constants[155])-exp(algebraic[148]/constants[155])) , True, (((constants[22]*algebraic[92])/constants[155])/(exp(algebraic[92]/constants[155])-1.00000))*(exp(algebraic[145]/constants[155])-exp(algebraic[148]/constants[155]))]) algebraic[157] = -algebraic[134]-algebraic[156] algebraic[158] = algebraic[156]+algebraic[134] algebraic[97] = constants[155]*log(constants[79]*states[37]) algebraic[99] = constants[155]*log(constants[85]*states[43]) algebraic[159] = constants[41]*(exp(algebraic[97]/constants[155])-exp((algebraic[99]+constants[146]*algebraic[80])/constants[155])) algebraic[101] = constants[155]*log(constants[80]*states[38]) algebraic[105] = constants[155]*log(constants[86]*states[44]) algebraic[160] = constants[43]*(exp(algebraic[101]/constants[155])-exp((algebraic[105]+constants[146]*algebraic[80])/constants[155])) algebraic[107] = constants[155]*log(constants[81]*states[39]) algebraic[109] = constants[155]*log(constants[87]*states[45]) algebraic[161] = constants[45]*(exp(algebraic[107]/constants[155])-exp((algebraic[109]+constants[146]*algebraic[80])/constants[155])) algebraic[112] = constants[155]*log(constants[82]*states[40]) algebraic[114] = constants[155]*log(constants[88]*states[46]) algebraic[162] = constants[42]*(exp(algebraic[112]/constants[155])-exp((algebraic[114]+constants[146]*algebraic[80])/constants[155])) algebraic[116] = constants[155]*log(constants[83]*states[41]) algebraic[163] = constants[44]*(exp(algebraic[116]/constants[155])-exp((algebraic[120]+constants[146]*algebraic[80])/constants[155])) algebraic[136] = constants[155]*log(constants[84]*states[42]) algebraic[164] = constants[46]*(exp(algebraic[136]/constants[155])-exp((algebraic[140]+constants[146]*algebraic[80])/constants[155])) algebraic[165] = (constants[130]*(algebraic[162]+algebraic[163]+algebraic[161]+algebraic[160]+algebraic[159]+algebraic[164])-algebraic[132])-algebraic[150] algebraic[172] = constants[29]*(exp((algebraic[97]+constants[148]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[101]+constants[152]*constants[96]*algebraic[76])/constants[155])) algebraic[173] = constants[33]*(exp((algebraic[97]+constants[149]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[107]+constants[153]*constants[96]*algebraic[76])/constants[155])) algebraic[166] = constants[23]*(exp((algebraic[97]+constants[147]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[112]+constants[151]*constants[96]*algebraic[76])/constants[155])) algebraic[175] = ((-algebraic[166]-algebraic[172])-algebraic[173])-algebraic[159] algebraic[174] = constants[37]*(exp((algebraic[101]+constants[150]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[107]+constants[154]*constants[96]*algebraic[76])/constants[155])) algebraic[168] = constants[25]*(exp((algebraic[107]+constants[147]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[136]+constants[151]*constants[96]*algebraic[76])/constants[155])) algebraic[177] = (-algebraic[168]+algebraic[174]+algebraic[173])-algebraic[161] algebraic[176] = constants[30]*(exp((algebraic[112]+constants[148]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[116]+constants[152]*constants[96]*algebraic[76])/constants[155])) algebraic[167] = constants[24]*(exp((algebraic[101]+constants[147]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[116]+constants[151]*constants[96]*algebraic[76])/constants[155])) algebraic[179] = ((-algebraic[167]-algebraic[174])+algebraic[176])-algebraic[160] algebraic[178] = constants[34]*(exp((algebraic[112]+constants[149]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[136]+constants[153]*constants[96]*algebraic[76])/constants[155])) algebraic[181] = ((-algebraic[178]+algebraic[166])-algebraic[176])-algebraic[162] algebraic[180] = constants[38]*(exp((algebraic[116]+constants[150]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[136]+constants[154]*constants[96]*algebraic[76])/constants[155])) algebraic[183] = ((algebraic[176]+algebraic[167])-algebraic[163])-algebraic[180] algebraic[184] = ((algebraic[178]+algebraic[168])-algebraic[164])+algebraic[180] algebraic[169] = constants[26]*(exp((algebraic[99]+constants[147]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[114]+constants[151]*constants[96]*algebraic[76])/constants[155])) algebraic[182] = constants[31]*(exp((algebraic[99]+constants[148]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[105]+constants[152]*constants[96]*algebraic[76])/constants[155])) algebraic[185] = constants[35]*(exp((algebraic[99]+constants[149]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[109]+constants[153]*constants[96]*algebraic[76])/constants[155])) algebraic[187] = ((-algebraic[169]-algebraic[182])-algebraic[185])+algebraic[159] algebraic[171] = constants[28]*(exp((algebraic[109]+constants[147]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[140]+constants[151]*constants[96]*algebraic[76])/constants[155])) algebraic[186] = constants[39]*(exp((algebraic[105]+constants[150]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[109]+constants[154]*constants[96]*algebraic[76])/constants[155])) algebraic[189] = -algebraic[171]+algebraic[186]+algebraic[185]+algebraic[161] algebraic[170] = constants[27]*(exp((algebraic[105]+constants[147]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[120]+constants[151]*constants[96]*algebraic[76])/constants[155])) algebraic[188] = constants[32]*(exp((algebraic[114]+constants[148]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[120]+constants[152]*constants[96]*algebraic[76])/constants[155])) algebraic[191] = (-algebraic[170]-algebraic[186])+algebraic[188]+algebraic[160] algebraic[190] = constants[36]*(exp((algebraic[114]+constants[149]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[140]+constants[153]*constants[96]*algebraic[76])/constants[155])) algebraic[193] = ((-algebraic[190]+algebraic[169])-algebraic[188])+algebraic[162] algebraic[192] = constants[40]*(exp((algebraic[120]+constants[150]*constants[96]*algebraic[76])/constants[155])-exp((algebraic[140]+constants[154]*constants[96]*algebraic[76])/constants[155])) algebraic[195] = (algebraic[188]+algebraic[170]+algebraic[163])-algebraic[192] algebraic[196] = algebraic[190]+algebraic[171]+algebraic[164]+algebraic[192] algebraic[73] = constants[130]*algebraic[70]*constants[96] algebraic[194] = algebraic[185]*(constants[153]-constants[149])+algebraic[186]*(constants[154]-constants[150])+algebraic[188]*(constants[152]-constants[148])+algebraic[190]*(constants[153]-constants[149])+algebraic[192]*(constants[154]-constants[150])+algebraic[182]*(constants[152]-constants[148])+algebraic[171]*(constants[151]-constants[147])+algebraic[170]*(constants[151]-constants[147])+algebraic[169]*(constants[151]-constants[147])+algebraic[180]*(constants[154]-constants[150])+algebraic[178]*(constants[153]-constants[149])+algebraic[176]*(constants[152]-constants[148])+algebraic[174]*(constants[154]-constants[150])+algebraic[173]*(constants[153]-constants[149])+algebraic[172]*(constants[152]-constants[148])+algebraic[168]*(constants[151]-constants[147])+algebraic[167]*(constants[151]-constants[147])+algebraic[166]*(constants[151]-constants[147]) algebraic[197] = constants[96]*((algebraic[194]-constants[130]*(algebraic[132]+algebraic[150]))-constants[131]*(algebraic[134]+algebraic[156])) algebraic[198] = (-algebraic[73]-algebraic[197])+algebraic[49] algebraic[2] = states[19]/constants[139] algebraic[3] = states[20]/constants[156] algebraic[4] = states[18]/constants[157] algebraic[5] = states[22]/constants[139] algebraic[6] = states[21]/constants[139] algebraic[7] = states[23]/constants[139] algebraic[77] = -algebraic[74] algebraic[102] = -constants[136]*algebraic[100] algebraic[118] = algebraic[115] algebraic[125] = algebraic[123] algebraic[137] = constants[137]*algebraic[135] algebraic[152] = ((constants[138]*algebraic[108]-algebraic[93])-constants[138]*algebraic[144])+algebraic[151] return algebraic def custom_piecewise(cases): """Compute result of a piecewise function""" return select(cases[0::2],cases[1::2]) def solve_model(): (legend_states, legend_algebraic, legend_voi, legend_constants) = createLegends() """Solve model with ODE solver""" from scipy.integrate import ode # Initialise constants and state variables. Multiply vertain enzyme concentrations with Gmult factor glabels = [ # 'q_Gs_init in component environment (fmol)', ] [igs, _, _] = find_indices(glabels, legend_constants, legend_states, legend_algebraic); Gmult = 1e0; # 5e1; (init_states, constants) = initConsts(Gmult, igs) # Set timespan to solve over voi = linspace(0, 1, 15000) # options for stimulus protocol stimLabels = ["pulse_start in component environment (second)", "pulse_end in component environment (second)", "pulseMag in component environment (fA)", "pulseHolding in component environment (fA)" ] sd = find_indices_from_labels(stimLabels, legend_constants) constants[sd['pulse_start']] = 1.1e-3 constants[sd['pulse_start']] = 3e-3 constants[sd['pulseMag']] = 1.888e7 constants[sd['pulseHolding']] = 0 # Construct ODE object to solve r = ode(computeRates) r.set_integrator('vode', method='bdf', atol=1e-07, rtol=1e-07, max_step=1e-6) 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) for glab in glabels: print('Gmult applied for ' + glab) return (voi, constants, states, algebraic, Gmult) def plot_model(voi, constants, states, algebraic, Gmult): """Plot variables against variable of integration""" (legend_states, legend_algebraic, legend_voi, legend_constants) = createLegends() labels = ['I_pulse in component environment (fA)', 'V_m in component environment (volt)', 'q_Cai in component environment (fmol)', 'q_Ca_SR in component environment (fmol)', 'q_Ca_di in component environment (fmol)', 'q_Cao in component environment (fmol)', 'q_Ca_TRPN in component environment (fmol)', 'q_TRPN in component environment (fmol)'] [_, i_st, i_alg] = find_indices(labels, legend_constants, legend_states, legend_algebraic) if i_st: plot_selected(i_st, voi, states, legend_voi, legend_states, ('Gmult=%g' % (Gmult)), int(ceil(sqrt(len(i_st))))) if i_alg: plot_selected(i_alg, voi, algebraic, legend_voi, legend_algebraic, ('Gmult=%g' % (Gmult)), int(ceil(sqrt(len(i_alg))))) # plot all q variables to confirm if we are SS if False: plt.figure() n = ceil(sqrt(size(states, 2))) for i in range(1, size(states, 2)): plt.subplot(n, n, i + 1) plt.plot(voi, states[:, i]) xlabel(legend_voi) title(legend_states[i]) plt.show() def find_indices(labels, legend_constants, legend_states, legend_algebraic): i_con = [] i_st = [] i_alg = [] for lab in labels: try: i_con.append(legend_constants.index(lab)) except: try: i_alg.append(legend_algebraic.index(lab)) except: i_st.append(legend_states.index(lab)) return (i_con, i_st, i_alg) def find_indices_from_labels(labels, legend_list): # return indices for a name in a dict, with the key being the name of the label, andits value as the index sd = dict() i_found = [legend_list.index(lab) for lab in labels] sd = {labels[i].split(' ')[0]: i_found[i] for i in range(len(labels))} return sd def plot_selected(ids, x, y, legend_x, legend_y, titlestr, ns): istart = 30; plt.figure(); for i, id in enumerate(ids): plt.subplot(ns, ns, i + 1) plt.plot(x[istart:-1], y[id, istart:-1]) plt.xlabel('time (s)') plt.legend([legend_y[id].split(' ')[0]]) # plt.tight_layout() plt.subplots_adjust(wspace=0.2, hspace=0.4) plt.suptitle(titlestr) if __name__ == "__main__": t = time.time() (voi, constants,states, algebraic,Gmult) = solve_model() print('elapsed = ', round(time.time() - t,3), ' s') plot_model(voi, constants,states, algebraic, Gmult)