# Size of variable arrays: sizeAlgebraic = 121 sizeStates = 31 sizeConstants = 134 from math import * from numpy import * def createLegends(): legend_states = [""] * sizeStates legend_rates = [""] * sizeStates legend_algebraic = [""] * sizeAlgebraic legend_voi = "" legend_constants = [""] * sizeConstants legend_voi = "time in component environment (second)" legend_constants[0] = "protocol in component environment (dimensionless)" legend_constants[1] = "L_totmax in component b1_AR_Gs_parameters (uM)" legend_constants[2] = "sum_b1_AR in component b1_AR_Gs_parameters (uM)" legend_constants[3] = "Gs_tot in component b1_AR_Gs_parameters (uM)" legend_constants[4] = "Kl in component b1_AR_Gs_parameters (uM)" legend_constants[5] = "Kr in component b1_AR_Gs_parameters (uM)" legend_constants[6] = "Kc in component b1_AR_Gs_parameters (uM)" legend_constants[7] = "k_bar_kp in component b1_AR_Gs_parameters (per_sec)" legend_constants[8] = "k_bar_km in component b1_AR_Gs_parameters (per_sec)" legend_constants[9] = "k_p_kap in component b1_AR_Gs_parameters (per_uM_per_sec)" legend_constants[10] = "k_p_kam in component b1_AR_Gs_parameters (per_sec)" legend_constants[11] = "k_g_act in component b1_AR_Gs_parameters (per_sec)" legend_constants[12] = "k_hyd in component b1_AR_Gs_parameters (per_sec)" legend_constants[13] = "k_reassoc in component b1_AR_Gs_parameters (per_uM_per_sec)" legend_constants[14] = "AC_tot in component cAMP_parameters (uM)" legend_constants[15] = "ATP in component cAMP_parameters (uM)" legend_constants[16] = "PDE_tot in component cAMP_parameters (uM)" legend_constants[17] = "IBMX_tot in component cAMP_parameters (uM)" legend_constants[18] = "Fsk_tot in component cAMP_parameters (uM)" legend_constants[19] = "k_ac_basal in component cAMP_parameters (per_sec)" legend_constants[20] = "k_ac_gsa in component cAMP_parameters (per_sec)" legend_constants[21] = "k_ac_fsk in component cAMP_parameters (per_sec)" legend_constants[22] = "k_pde in component cAMP_parameters (per_sec)" legend_constants[23] = "Km_basal in component cAMP_parameters (uM)" legend_constants[24] = "Km_gsa in component cAMP_parameters (uM)" legend_constants[25] = "Km_fsk in component cAMP_parameters (uM)" legend_constants[26] = "Km_pde in component cAMP_parameters (uM)" legend_constants[27] = "K_gsa in component cAMP_parameters (uM)" legend_constants[28] = "K_fsk in component cAMP_parameters (uM)" legend_constants[29] = "Ki_ibmx in component cAMP_parameters (uM)" legend_constants[30] = "PKAI_tot in component PKA_parameters (uM)" legend_constants[31] = "PKAII_tot in component PKA_parameters (uM)" legend_constants[32] = "PKI_tot in component PKA_parameters (uM)" legend_constants[33] = "K_a in component PKA_parameters (uM)" legend_constants[34] = "K_b in component PKA_parameters (uM)" legend_constants[35] = "K_d in component PKA_parameters (uM)" legend_constants[36] = "Ki_pki in component PKA_parameters (uM)" legend_constants[37] = "PLB_tot in component PLB_parameters (uM)" legend_constants[38] = "PP1_tot in component PLB_parameters (uM)" legend_constants[39] = "Inhib1_tot in component PLB_parameters (uM)" legend_constants[40] = "k_pka_plb in component PLB_parameters (per_sec)" legend_constants[41] = "Km_pka_plb in component PLB_parameters (uM)" legend_constants[42] = "k_pp1_plb in component PLB_parameters (per_sec)" legend_constants[43] = "Km_pp1_plb in component PLB_parameters (uM)" legend_constants[44] = "k_pka_i1 in component PLB_parameters (per_sec)" legend_constants[45] = "Km_pka_i1 in component PLB_parameters (uM)" legend_constants[46] = "Vmax_pp2a_i1 in component PLB_parameters (uM_per_sec)" legend_constants[47] = "Km_pp2a_i1 in component PLB_parameters (uM)" legend_constants[48] = "Ki_inhib1 in component PLB_parameters (uM)" legend_constants[49] = "epsilon in component LCC_parameters (dimensionless)" legend_constants[50] = "LCC_tot in component LCC_parameters (uM)" legend_constants[51] = "PP1_lcc_tot in component LCC_parameters (uM)" legend_constants[52] = "PP2A_lcc_tot in component LCC_parameters (uM)" legend_constants[53] = "k_pka_lcc in component LCC_parameters (per_sec)" legend_constants[54] = "Km_pka_lcc in component LCC_parameters (uM)" legend_constants[55] = "k_pp1_lcc in component LCC_parameters (per_sec)" legend_constants[56] = "Km_pp1_lcc in component LCC_parameters (uM)" legend_constants[57] = "k_pp2a_lcc in component LCC_parameters (per_sec)" legend_constants[58] = "Km_pp2a_lcc in component LCC_parameters (uM)" legend_constants[59] = "V_myo in component EC_Coupling_Parameters (uL)" legend_constants[60] = "V_nsr in component EC_Coupling_Parameters (uL)" legend_constants[61] = "V_jsr in component EC_Coupling_Parameters (uL)" legend_constants[62] = "A_Cap in component EC_Coupling_Parameters (cm2)" legend_constants[63] = "Temp in component EC_Coupling_Parameters (kelvin)" legend_constants[64] = "Na_ext in component EC_Coupling_Parameters (mM)" legend_constants[65] = "K_ext in component EC_Coupling_Parameters (mM)" legend_constants[66] = "Ca_ext in component EC_Coupling_Parameters (mM)" legend_constants[67] = "G_Na in component EC_Coupling_Parameters (mS_per_uF)" legend_constants[68] = "G_to in component EC_Coupling_Parameters (mS_per_uF)" legend_constants[69] = "G_ss in component EC_Coupling_Parameters (mS_per_uF)" legend_constants[70] = "G_Ki_bar in component EC_Coupling_Parameters (mS_per_uF)" legend_constants[71] = "G_Kp in component EC_Coupling_Parameters (mS_per_uF)" legend_constants[72] = "f in component EC_Coupling_Parameters (per_sec)" legend_constants[73] = "g in component EC_Coupling_Parameters (per_sec)" legend_constants[74] = "gamma_o in component EC_Coupling_Parameters (per_mM_per_sec)" legend_constants[75] = "omega in component EC_Coupling_Parameters (per_sec)" legend_constants[76] = "p_Ca in component EC_Coupling_Parameters (cm_per_sec)" legend_constants[77] = "p_K in component EC_Coupling_Parameters (cm_per_sec)" legend_constants[78] = "N_lcc in component EC_Coupling_Parameters (dimensionless)" legend_constants[79] = "I_Ca05 in component EC_Coupling_Parameters (uA_per_uF)" legend_constants[80] = "k_NaCa in component EC_Coupling_Parameters (uA_per_uF)" legend_constants[81] = "Km_Na in component EC_Coupling_Parameters (mM)" legend_constants[82] = "Km_Ca in component EC_Coupling_Parameters (mM)" legend_constants[83] = "k_sat in component EC_Coupling_Parameters (dimensionless)" legend_constants[84] = "eta in component EC_Coupling_Parameters (dimensionless)" legend_constants[85] = "I_bar_NaK in component EC_Coupling_Parameters (uA_per_uF)" legend_constants[86] = "Km_Nai in component EC_Coupling_Parameters (mM)" legend_constants[87] = "Km_Ko in component EC_Coupling_Parameters (mM)" legend_constants[88] = "I_bar_PCa in component EC_Coupling_Parameters (uA_per_uF)" legend_constants[89] = "Km_PCa in component EC_Coupling_Parameters (mM)" legend_constants[90] = "G_CaB in component EC_Coupling_Parameters (uA_per_uF)" legend_constants[91] = "G_NaB in component EC_Coupling_Parameters (uA_per_uF)" legend_constants[92] = "Pns in component EC_Coupling_Parameters (dimensionless)" legend_constants[93] = "Km_NS in component EC_Coupling_Parameters (mM)" legend_constants[94] = "I_up_bar in component EC_Coupling_Parameters (mM_per_sec)" legend_constants[95] = "Km_up0 in component EC_Coupling_Parameters (mM)" legend_constants[96] = "NSR_bar in component EC_Coupling_Parameters (mM)" legend_constants[97] = "tau_on in component EC_Coupling_Parameters (second)" legend_constants[98] = "tau_off in component EC_Coupling_Parameters (second)" legend_constants[99] = "G_max_rel in component EC_Coupling_Parameters (mM_per_sec)" legend_constants[100] = "d_Cai_th in component EC_Coupling_Parameters (mM)" legend_constants[101] = "Km_rel in component EC_Coupling_Parameters (mM)" legend_constants[102] = "CSQN_th in component EC_Coupling_Parameters (mM)" legend_constants[103] = "CSQN_bar in component EC_Coupling_Parameters (mM)" legend_constants[104] = "Km_CSQN in component EC_Coupling_Parameters (mM)" legend_constants[105] = "tau_tr in component EC_Coupling_Parameters (second)" legend_constants[106] = "TRPN_bar in component EC_Coupling_Parameters (mM)" legend_constants[107] = "CMDN_bar in component EC_Coupling_Parameters (mM)" legend_constants[108] = "INDO_bar in component EC_Coupling_Parameters (mM)" legend_constants[109] = "Km_TRPN in component EC_Coupling_Parameters (mM)" legend_constants[110] = "Km_CMDN in component EC_Coupling_Parameters (mM)" legend_constants[111] = "Km_INDO in component EC_Coupling_Parameters (mM)" legend_algebraic[61] = "LR in component b1_AR_module (uM)" legend_algebraic[62] = "LRG in component b1_AR_module (uM)" legend_algebraic[63] = "RG in component b1_AR_module (uM)" legend_algebraic[79] = "BARK_DESENS in component b1_AR_module (uM_per_sec)" legend_algebraic[11] = "BARK_RESENS in component b1_AR_module (uM_per_sec)" legend_algebraic[87] = "PKA_DESENS in component b1_AR_module (uM_per_sec)" legend_algebraic[22] = "PKA_RESENS in component b1_AR_module (uM_per_sec)" legend_algebraic[80] = "G_ACT in component b1_AR_module (uM_per_sec)" legend_algebraic[24] = "HYD in component b1_AR_module (uM_per_sec)" legend_algebraic[26] = "REASSOC in component b1_AR_module (uM_per_sec)" legend_algebraic[64] = "L in component b1_AR_module (uM)" legend_algebraic[65] = "R in component b1_AR_module (uM)" legend_algebraic[66] = "Gs in component b1_AR_module (uM)" legend_states[0] = "b1_AR_tot in component b1_AR_module (uM)" legend_states[1] = "b1_AR_d in component b1_AR_module (uM)" legend_states[2] = "b1_AR_p in component b1_AR_module (uM)" legend_states[3] = "Gs_agtp_tot in component b1_AR_module (uM)" legend_states[4] = "Gs_agdp in component b1_AR_module (uM)" legend_states[5] = "Gs_bg in component b1_AR_module (uM)" legend_algebraic[67] = "PKAC_I in component PKA_module (uM)" legend_algebraic[68] = "cAMP in component PKA_module (uM)" legend_algebraic[41] = "Gsa_GTP in component cAMP_module (uM)" legend_algebraic[49] = "Fsk in component cAMP_module (uM)" legend_algebraic[42] = "AC in component cAMP_module (uM)" legend_constants[131] = "PDE in component cAMP_module (uM)" legend_constants[132] = "IBMX in component cAMP_module (uM)" legend_states[6] = "cAMP_tot in component cAMP_module (uM)" legend_algebraic[43] = "Gsa_GTP_AC in component cAMP_module (uM)" legend_algebraic[50] = "Fsk_AC in component cAMP_module (uM)" legend_algebraic[47] = "AC_ACT_GSA in component cAMP_module (uM_per_sec)" legend_algebraic[45] = "AC_ACT_BASAL in component cAMP_module (uM_per_sec)" legend_algebraic[52] = "AC_ACT_FSK in component cAMP_module (uM_per_sec)" legend_algebraic[81] = "PDE_ACT in component cAMP_module (uM_per_sec)" legend_constants[133] = "PDE_IBMX in component cAMP_module (uM)" legend_algebraic[69] = "PKI in component PKA_module (uM)" legend_algebraic[70] = "A2RC_I in component PKA_module (uM)" legend_algebraic[71] = "A2R_I in component PKA_module (uM)" legend_algebraic[72] = "A2RC_II in component PKA_module (uM)" legend_algebraic[73] = "A2R_II in component PKA_module (uM)" legend_algebraic[74] = "ARC_I in component PKA_module (uM)" legend_algebraic[75] = "ARC_II in component PKA_module (uM)" legend_algebraic[76] = "PKA_temp in component PKA_module (uM)" legend_algebraic[77] = "PKAC_II in component PKA_module (uM)" legend_constants[112] = "zed in component PKA_module (dimensionless)" legend_constants[113] = "zed1 in component PKA_module (dimensionless)" legend_algebraic[28] = "PLB in component PLB_module (uM)" legend_algebraic[82] = "PLB_PHOSPH in component PLB_module (uM_per_sec)" legend_algebraic[59] = "PLB_DEPHOSPH in component PLB_module (uM_per_sec)" legend_algebraic[29] = "Inhib1 in component PLB_module (uM)" legend_algebraic[54] = "Inhib1p_PP1 in component PLB_module (uM)" legend_algebraic[83] = "Inhib1_PHOSPH in component PLB_module (uM_per_sec)" legend_algebraic[31] = "Inhib1_DEPHOSPH in component PLB_module (uM_per_sec)" legend_states[7] = "PLB_p in component PLB_module (uM)" legend_states[8] = "Inhib1p_tot in component PLB_module (uM)" legend_algebraic[55] = "Inhib1p in component PLB_module (uM)" legend_algebraic[56] = "PP1 in component PLB_module (uM)" legend_algebraic[0] = "frac_PLB_p in component PLB_module (dimensionless)" legend_algebraic[30] = "frac_PLB in component PLB_module (dimensionless)" legend_constants[124] = "frac_PLB_o in component PLB_module (dimensionless)" legend_algebraic[33] = "LCCa in component LCC_module (uM)" legend_algebraic[84] = "LCCa_PHOSPH in component LCC_module (uM_per_sec)" legend_algebraic[35] = "LCCa_DEPHOSPH in component LCC_module (uM_per_sec)" legend_algebraic[37] = "LCCb in component LCC_module (uM)" legend_algebraic[85] = "LCCb_PHOSPH in component LCC_module (uM_per_sec)" legend_algebraic[39] = "LCCb_DEPHOSPH in component LCC_module (uM_per_sec)" legend_states[9] = "LCCa_p in component LCC_module (uM)" legend_states[10] = "LCCb_p in component LCC_module (uM)" legend_algebraic[1] = "frac_LCCa_p in component LCC_module (dimensionless)" legend_constants[125] = "frac_LCCa_po in component LCC_module (dimensionless)" legend_algebraic[32] = "frac_LCCb_p in component LCC_module (dimensionless)" legend_constants[126] = "frac_LCCb_po in component LCC_module (dimensionless)" legend_algebraic[34] = "E_Na in component Nernst_Potentials (mV)" legend_algebraic[36] = "E_K in component Nernst_Potentials (mV)" legend_algebraic[38] = "E_Ca in component Nernst_Potentials (mV)" legend_constants[127] = "E_Cl in component Nernst_Potentials (mV)" legend_constants[114] = "R in component Nernst_Potentials (joules_per_kmole_kelvin)" legend_constants[115] = "Frdy in component Nernst_Potentials (coulombs_per_mole)" legend_constants[128] = "FoRT in component Nernst_Potentials (per_mV)" legend_constants[116] = "z_Na in component Nernst_Potentials (dimensionless)" legend_constants[117] = "z_K in component Nernst_Potentials (dimensionless)" legend_constants[118] = "z_Ca in component Nernst_Potentials (dimensionless)" legend_states[11] = "Na_i in component Ion_Concentrations_and_Membrane_Potential (mM)" legend_states[12] = "K_i in component Ion_Concentrations_and_Membrane_Potential (mM)" legend_states[13] = "Ca_i in component Ion_Concentrations_and_Membrane_Potential (mM)" legend_algebraic[2] = "am in component Fast_Na_Current (per_sec)" legend_algebraic[12] = "bm in component Fast_Na_Current (per_sec)" legend_algebraic[3] = "ah in component Fast_Na_Current (per_sec)" legend_algebraic[4] = "aj in component Fast_Na_Current (per_sec)" legend_algebraic[13] = "bh in component Fast_Na_Current (per_sec)" legend_algebraic[14] = "bj in component Fast_Na_Current (per_sec)" legend_states[14] = "m in component Fast_Na_Current (dimensionless)" legend_states[15] = "h in component Fast_Na_Current (dimensionless)" legend_states[16] = "j in component Fast_Na_Current (dimensionless)" legend_algebraic[40] = "I_Na in component Fast_Na_Current (uA_per_uF)" legend_states[17] = "V_m in component Ion_Concentrations_and_Membrane_Potential (mV)" legend_algebraic[5] = "a_lcc in component L_Type_Calcium_Current (per_sec)" legend_algebraic[15] = "b_lcc in component L_Type_Calcium_Current (per_sec)" legend_algebraic[17] = "f_lcc in component L_Type_Calcium_Current (per_sec)" legend_algebraic[6] = "y_lcc_inf in component L_Type_Calcium_Current (dimensionless)" legend_algebraic[16] = "tau_y_lcc in component L_Type_Calcium_Current (second)" legend_algebraic[23] = "gamma in component L_Type_Calcium_Current (per_sec)" legend_algebraic[25] = "v_gamma in component L_Type_Calcium_Current (per_sec)" legend_algebraic[27] = "v_omega in component L_Type_Calcium_Current (per_sec)" legend_states[18] = "v in component L_Type_Calcium_Current (dimensionless)" legend_states[19] = "w in component L_Type_Calcium_Current (dimensionless)" legend_states[20] = "x in component L_Type_Calcium_Current (dimensionless)" legend_states[21] = "y in component L_Type_Calcium_Current (dimensionless)" legend_states[22] = "z in component L_Type_Calcium_Current (dimensionless)" legend_algebraic[44] = "i_bar_Ca in component L_Type_Calcium_Current (uA_per_uF)" legend_algebraic[46] = "i_bar_K in component L_Type_Calcium_Current (uA_per_uF)" legend_algebraic[48] = "f_avail in component L_Type_Calcium_Current (dimensionless)" legend_algebraic[51] = "I_Ca in component L_Type_Calcium_Current (uA_per_uF)" legend_algebraic[53] = "I_CaK in component L_Type_Calcium_Current (uA_per_uF)" legend_algebraic[57] = "I_Ca_tot in component L_Type_Calcium_Current (uA_per_uF)" legend_constants[119] = "alpha in component L_Type_Calcium_Current (ibar_to_I)" legend_algebraic[7] = "r_toss in component Transient_Outward_K_Current (dimensionless)" legend_algebraic[8] = "s_toss in component Transient_Outward_K_Current (dimensionless)" legend_algebraic[18] = "tau_r_to in component Transient_Outward_K_Current (second)" legend_algebraic[19] = "tau_s_to in component Transient_Outward_K_Current (second)" legend_algebraic[20] = "tau_ss_to in component Transient_Outward_K_Current (second)" legend_states[23] = "r_to in component Transient_Outward_K_Current (dimensionless)" legend_states[24] = "s_to in component Transient_Outward_K_Current (dimensionless)" legend_states[25] = "ss_to in component Transient_Outward_K_Current (dimensionless)" legend_algebraic[58] = "I_to in component Transient_Outward_K_Current (uA_per_uF)" legend_algebraic[9] = "r_ss_inf in component Steady_State_K_Current (dimensionless)" legend_algebraic[21] = "tau_r_ss in component Steady_State_K_Current (second)" legend_algebraic[10] = "s_ss_inf in component Steady_State_K_Current (dimensionless)" legend_constants[129] = "tau_s_ss in component Steady_State_K_Current (second)" legend_states[26] = "r_ss in component Steady_State_K_Current (dimensionless)" legend_states[27] = "s_ss in component Steady_State_K_Current (dimensionless)" legend_algebraic[60] = "I_ss in component Steady_State_K_Current (uA_per_uF)" legend_algebraic[78] = "a_Ki in component Time_Independent_K_Current (dimensionless)" legend_algebraic[86] = "b_Ki in component Time_Independent_K_Current (dimensionless)" legend_algebraic[88] = "Ki_ss in component Time_Independent_K_Current (dimensionless)" legend_algebraic[89] = "I_Ki in component Time_Independent_K_Current (uA_per_uF)" legend_algebraic[90] = "Kp in component Plateau_K_Current (dimensionless)" legend_algebraic[91] = "I_Kp in component Plateau_K_Current (uA_per_uF)" legend_algebraic[92] = "s4 in component Na_Ca_Exchanger_Current (dimensionless)" legend_algebraic[93] = "s5 in component Na_Ca_Exchanger_Current (dimensionless)" legend_algebraic[94] = "I_NCX in component Na_Ca_Exchanger_Current (uA_per_uF)" legend_constants[130] = "sigma in component Na_K_Pump_Current (dimensionless)" legend_algebraic[95] = "f_NaK in component Na_K_Pump_Current (dimensionless)" legend_algebraic[96] = "I_NaK in component Na_K_Pump_Current (uA_per_uF)" legend_algebraic[97] = "I_PCa in component Sarcolemmal_Ca_Pump_Current (uA_per_uF)" legend_algebraic[98] = "I_CaB in component Ca_Background_Current (uA_per_uF)" legend_algebraic[99] = "I_NaB in component Na_Background_Current (uA_per_uF)" legend_algebraic[100] = "I_Na_tot in component Total_Membrane_Currents (uA_per_uF)" legend_algebraic[101] = "I_K_tot in component Total_Membrane_Currents (uA_per_uF)" legend_algebraic[102] = "I_Ca_tot in component Total_Membrane_Currents (uA_per_uF)" legend_algebraic[103] = "t_rel_sub in component Calcium_Induced_Calcium_Release (second)" legend_algebraic[105] = "ryr_on in component Calcium_Induced_Calcium_Release (dimensionless)" legend_algebraic[107] = "ryr_off in component Calcium_Induced_Calcium_Release (dimensionless)" legend_algebraic[109] = "g_rel in component Calcium_Induced_Calcium_Release (per_sec)" legend_algebraic[110] = "I_rel in component Calcium_Induced_Calcium_Release (mM_per_sec)" legend_states[28] = "Ca_jsr in component Other_SR_Fluxes_and_Concentrations (mM)" legend_states[29] = "trel in component Ion_Concentrations_and_Membrane_Potential (second)" legend_algebraic[111] = "Km_up in component Other_SR_Fluxes_and_Concentrations (mM)" legend_algebraic[112] = "I_up in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)" legend_algebraic[113] = "I_leak in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)" legend_algebraic[114] = "I_tr in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)" legend_algebraic[116] = "B_jsr in component Other_SR_Fluxes_and_Concentrations (dimensionless)" legend_states[30] = "Ca_nsr in component Other_SR_Fluxes_and_Concentrations (mM)" legend_algebraic[118] = "SR_content in component Other_SR_Fluxes_and_Concentrations (uM)" legend_algebraic[115] = "b_trpn in component Cytoplasmic_Calcium_Buffering (dimensionless)" legend_algebraic[117] = "b_cmdn in component Cytoplasmic_Calcium_Buffering (dimensionless)" legend_algebraic[119] = "b_indo in component Cytoplasmic_Calcium_Buffering (dimensionless)" legend_algebraic[120] = "B_myo in component Cytoplasmic_Calcium_Buffering (dimensionless)" legend_algebraic[108] = "I_app in component Ion_Concentrations_and_Membrane_Potential (uA_per_uF)" legend_constants[120] = "V_hold in component Ion_Concentrations_and_Membrane_Potential (mV)" legend_constants[121] = "V_test in component Ion_Concentrations_and_Membrane_Potential (mV)" legend_algebraic[104] = "V_clamp in component Ion_Concentrations_and_Membrane_Potential (mV)" legend_constants[122] = "R_clamp in component Ion_Concentrations_and_Membrane_Potential (ms)" legend_constants[123] = "lambda in component Ion_Concentrations_and_Membrane_Potential (nFC_per_sAcm2)" legend_algebraic[106] = "I_pace in component Ion_Concentrations_and_Membrane_Potential (uA_per_uF)" legend_rates[0] = "d/dt b1_AR_tot in component b1_AR_module (uM)" legend_rates[1] = "d/dt b1_AR_d in component b1_AR_module (uM)" legend_rates[2] = "d/dt b1_AR_p in component b1_AR_module (uM)" legend_rates[3] = "d/dt Gs_agtp_tot in component b1_AR_module (uM)" legend_rates[4] = "d/dt Gs_agdp in component b1_AR_module (uM)" legend_rates[5] = "d/dt Gs_bg in component b1_AR_module (uM)" legend_rates[6] = "d/dt cAMP_tot in component cAMP_module (uM)" legend_rates[7] = "d/dt PLB_p in component PLB_module (uM)" legend_rates[8] = "d/dt Inhib1p_tot in component PLB_module (uM)" legend_rates[9] = "d/dt LCCa_p in component LCC_module (uM)" legend_rates[10] = "d/dt LCCb_p in component LCC_module (uM)" legend_rates[14] = "d/dt m in component Fast_Na_Current (dimensionless)" legend_rates[15] = "d/dt h in component Fast_Na_Current (dimensionless)" legend_rates[16] = "d/dt j in component Fast_Na_Current (dimensionless)" legend_rates[18] = "d/dt v in component L_Type_Calcium_Current (dimensionless)" legend_rates[19] = "d/dt w in component L_Type_Calcium_Current (dimensionless)" legend_rates[20] = "d/dt x in component L_Type_Calcium_Current (dimensionless)" legend_rates[21] = "d/dt y in component L_Type_Calcium_Current (dimensionless)" legend_rates[22] = "d/dt z in component L_Type_Calcium_Current (dimensionless)" legend_rates[23] = "d/dt r_to in component Transient_Outward_K_Current (dimensionless)" legend_rates[24] = "d/dt s_to in component Transient_Outward_K_Current (dimensionless)" legend_rates[25] = "d/dt ss_to in component Transient_Outward_K_Current (dimensionless)" legend_rates[26] = "d/dt r_ss in component Steady_State_K_Current (dimensionless)" legend_rates[27] = "d/dt s_ss in component Steady_State_K_Current (dimensionless)" legend_rates[30] = "d/dt Ca_nsr in component Other_SR_Fluxes_and_Concentrations (mM)" legend_rates[28] = "d/dt Ca_jsr in component Other_SR_Fluxes_and_Concentrations (mM)" legend_rates[11] = "d/dt Na_i in component Ion_Concentrations_and_Membrane_Potential (mM)" legend_rates[12] = "d/dt K_i in component Ion_Concentrations_and_Membrane_Potential (mM)" legend_rates[13] = "d/dt Ca_i in component Ion_Concentrations_and_Membrane_Potential (mM)" legend_rates[17] = "d/dt V_m in component Ion_Concentrations_and_Membrane_Potential (mV)" legend_rates[29] = "d/dt trel in component Ion_Concentrations_and_Membrane_Potential (second)" return (legend_states, legend_algebraic, legend_voi, legend_constants) def initConsts(): constants = [0.0] * sizeConstants; states = [0.0] * sizeStates; constants[0] = 1 constants[1] = 0.01 constants[2] = 0.0132 constants[3] = 3.83 constants[4] = 0.285 constants[5] = 0.062 constants[6] = 33 constants[7] = 1.1e-3 constants[8] = 2.2e-3 constants[9] = 3.6e-3 constants[10] = 2.2e-3 constants[11] = 16 constants[12] = 0.8 constants[13] = 1.21e3 constants[14] = 49.7e-3 constants[15] = 5e3 constants[16] = 38.9e-3 constants[17] = 0 constants[18] = 0 constants[19] = 0.2 constants[20] = 8.5 constants[21] = 7.3 constants[22] = 5 constants[23] = 1.03e3 constants[24] = 315 constants[25] = 860 constants[26] = 1.3 constants[27] = 0.4 constants[28] = 44 constants[29] = 30 constants[30] = 0.59 constants[31] = 0.025 constants[32] = 0.18 constants[33] = 9.14 constants[34] = 1.64 constants[35] = 4.375 constants[36] = 0.2e-3 constants[37] = 106 constants[38] = 0.89 constants[39] = 0.3 constants[40] = 54 constants[41] = 21 constants[42] = 8.5 constants[43] = 7 constants[44] = 60 constants[45] = 1 constants[46] = 14 constants[47] = 1 constants[48] = 1e-3 constants[49] = 10 constants[50] = 25e-3 constants[51] = 25e-3 constants[52] = 25e-3 constants[53] = 54 constants[54] = 21 constants[55] = 8.52 constants[56] = 3 constants[57] = 10.1 constants[58] = 3 constants[59] = 20.8e-6 constants[60] = 9.88e-7 constants[61] = 9.3e-8 constants[62] = 1.534e-4 constants[63] = 310 constants[64] = 140 constants[65] = 5.4 constants[66] = 1.8 constants[67] = 8 constants[68] = 0.35 constants[69] = 0.07 constants[70] = 0.24 constants[71] = 0.008 constants[72] = 300 constants[73] = 2e3 constants[74] = 5187.5 constants[75] = 10 constants[76] = 1.7469e-8 constants[77] = 3.234e-11 constants[78] = 3e5 constants[79] = -0.458 constants[80] = 1483 constants[81] = 87.5 constants[82] = 1.38 constants[83] = 0.1 constants[84] = 0.35 constants[85] = 1.1 constants[86] = 10 constants[87] = 1.5 constants[88] = 1.15 constants[89] = 0.5e-3 constants[90] = 2.8e-3 constants[91] = 1.18e-3 constants[92] = 0 constants[93] = 1.2e-3 constants[94] = 4.7 constants[95] = 3e-4 constants[96] = 15 constants[97] = 2e-3 constants[98] = 2e-3 constants[99] = 60e3 constants[100] = 0.18e-3 constants[101] = 0.8e-3 constants[102] = 8.75 constants[103] = 15 constants[104] = 0.8 constants[105] = 5.7e-4 constants[106] = 0.07 constants[107] = 0.05 constants[108] = 0.07 constants[109] = 0.5128e-3 constants[110] = 2.38e-3 constants[111] = 8.44e-4 states[0] = 0.01205 states[1] = 0 states[2] = 1.154e-3 states[3] = 0.02505 states[4] = 6.446e-4 states[5] = 0.02569 states[6] = 0.8453 constants[112] = 0 constants[113] = 0 states[7] = 4.105 states[8] = 0.0526 states[9] = 5.103e-3 states[10] = 5.841e-3 constants[114] = 8314 constants[115] = 96485 constants[116] = 1 constants[117] = 1 constants[118] = 2 states[11] = 16 states[12] = 145 states[13] = 1.58e-4 states[14] = 1.4e-3 states[15] = 0.99 states[16] = 0.99 states[17] = -85.66 states[18] = 0 states[19] = 0 states[20] = 0.13 states[21] = 0.96 states[22] = 0.92 constants[119] = 10e5 states[23] = 1.4e-3 states[24] = 1 states[25] = 0.613 states[26] = 198e-3 states[27] = 0.43 states[28] = 1.92 states[29] = 0.9 states[30] = 1.92 constants[120] = -40 constants[121] = -10 constants[122] = 0.02 constants[123] = -1e3 constants[124] = 0.961300 constants[125] = 0.204100 constants[126] = 0.233600 constants[127] = -40.0000 constants[128] = (constants[115]/constants[114])/constants[63] constants[129] = 2.10000 constants[130] = (exp(constants[64]/67.3000)-1.00000)/7.00000 rootfind_0(voi, constants, rates, states, algebraic) return (states, constants) def computeRates(voi, states, constants): rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic algebraic[10] = 1.00000/(1.00000+exp((states[17]+87.5000)/10.3000)) rates[27] = (algebraic[10]-states[27])/constants[129] algebraic[2] = (0.320000*(states[17]+47.1300))/(1.00000-exp(-0.100000*(states[17]+47.1300))) algebraic[12] = 0.0800000*exp(-states[17]/11.0000) rates[14] = 1000.00*(algebraic[2]*(1.00000-states[14])-algebraic[12]*states[14]) algebraic[3] = custom_piecewise([greater_equal(states[17] , -40.0000), 0.00000 , True, 0.135000*exp((80.0000+states[17])/-6.80000)]) algebraic[13] = custom_piecewise([greater_equal(states[17] , -40.0000), 1.00000/(0.130000*(1.00000+exp(-(states[17]+10.6600)/11.1000))) , True, 3.56000*exp(0.0790000*states[17])+310000.*exp(0.350000*states[17])]) rates[15] = 1000.00*(algebraic[3]*(1.00000-states[15])-algebraic[13]*states[15]) algebraic[4] = custom_piecewise([greater_equal(states[17] , -40.0000), 0.00000 , True, ((-127140.*exp(0.244400*states[17])-3.47400e-05*exp(-0.0439100*states[17]))*(states[17]+37.7800))/(1.00000+exp(0.311000*(states[17]+79.2300)))]) algebraic[14] = custom_piecewise([greater_equal(states[17] , -40.0000), (0.300000*exp(-2.57500e-07*states[17]))/(1.00000+exp(-0.100000*(states[17]+32.0000))) , True, (0.121200*exp(-0.0105200*states[17]))/(1.00000+exp(-0.137800*(states[17]+40.1400)))]) rates[16] = 1000.00*(algebraic[4]*(1.00000-states[16])-algebraic[14]*states[16]) algebraic[5] = 400.000*exp((states[17]+2.00000)/10.0000) algebraic[15] = 50.0000*exp((-1.00000*(states[17]+2.00000))/13.0000) rates[18] = algebraic[5]*(1.00000-states[18])-algebraic[15]*states[18] rates[19] = 2.00000*algebraic[5]*(1.00000-states[19])-(algebraic[15]/2.00000)*states[19] algebraic[1] = states[9]/constants[50] algebraic[17] = constants[72]*((0.375000*algebraic[1])/constants[125]+0.625000) rates[20] = algebraic[17]*(1.00000-states[20])-constants[73]*states[20] algebraic[6] = 1.00000/(1.00000+exp((states[17]+55.0000)/7.50000))+0.100000/(1.00000+exp((-states[17]+21.0000)/6.00000)) algebraic[16] = 0.0200000+0.300000/(1.00000+exp((states[17]+30.0000)/9.50000)) rates[21] = (algebraic[6]-states[21])/algebraic[16] algebraic[7] = 1.00000/(1.00000+exp((states[17]+10.6000)/-11.4200)) algebraic[18] = 1.00000/(45.1600*exp(0.0357700*(states[17]+50.0000))+98.9000*exp(-0.100000*(states[17]+38.0000))) rates[23] = (algebraic[7]-states[23])/algebraic[18] algebraic[8] = 1.00000/(1.00000+exp((states[17]+43.5000)/6.88410)) algebraic[19] = 0.350000*exp(-1.00000*(power((states[17]+70.0000)/15.0000, 2.00000)))+0.0350000 rates[24] = (algebraic[8]-states[24])/algebraic[19] algebraic[20] = 3.70000*exp(-1.00000*(power((states[17]+70.0000)/30.0000, 2.00000)))+0.0350000 rates[25] = (algebraic[8]-states[25])/algebraic[20] algebraic[9] = 1.00000/(1.00000+exp((states[17]+11.5000)/-11.8200)) algebraic[21] = 10.0000/(45.1600*exp(0.0357700*(states[17]+50.0000))+98.9000*exp(-0.100000*(states[17]+38.0000))) rates[26] = (algebraic[9]-states[26])/algebraic[21] algebraic[24] = constants[12]*states[3] algebraic[26] = constants[13]*states[4]*states[5] rates[4] = algebraic[24]-algebraic[26] algebraic[23] = constants[74]*states[13] algebraic[25] = algebraic[23]*(power(1.00000-states[18], 4.00000)+2.00000*states[18]*(power(1.00000-states[18], 3.00000))+4.00000*(power(states[18], 2.00000))*(power(1.00000-states[18], 2.00000))+8.00000*(power(states[18], 3.00000))*(1.00000-states[18])+16.0000*(power(states[18], 4.00000))*(1.00000-algebraic[17]/constants[73])) algebraic[27] = constants[75]*(power(1.00000-states[19], 4.00000)+0.500000*states[19]*(power(1.00000-states[19], 3.00000))+0.250000*(power(states[19], 2.00000))*(power(1.00000-states[19], 2.00000))+0.125000*(power(states[19], 3.00000))*(1.00000-states[19])+0.0625000*(power(states[19], 4.00000))) rates[22] = algebraic[27]*(1.00000-states[22])-algebraic[25]*states[22] rootfind_1(voi, constants, rates, states, algebraic) algebraic[79] = constants[7]*(algebraic[61]+algebraic[62]) algebraic[11] = constants[8]*states[1] rates[1] = algebraic[79]-algebraic[11] algebraic[80] = constants[11]*(algebraic[63]+algebraic[62]) rates[3] = algebraic[80]-algebraic[24] rates[5] = algebraic[80]-algebraic[26] rootfind_2(voi, constants, rates, states, algebraic) algebraic[47] = (constants[20]*algebraic[43]*constants[15])/(constants[24]+constants[15]) algebraic[45] = (constants[19]*algebraic[42]*constants[15])/(constants[23]+constants[15]) rootfind_3(voi, constants, rates, states, algebraic) algebraic[52] = (constants[21]*algebraic[50]*constants[15])/(constants[25]+constants[15]) algebraic[81] = (constants[22]*constants[131]*algebraic[68])/(constants[26]+algebraic[68]) rates[6] = (algebraic[45]+algebraic[47]+algebraic[52])-algebraic[81] algebraic[28] = constants[37]-states[7] algebraic[82] = (constants[40]*algebraic[67]*algebraic[28])/(constants[41]+algebraic[28]) rootfind_4(voi, constants, rates, states, algebraic) algebraic[59] = (constants[42]*algebraic[56]*states[7])/(constants[43]+states[7]) rates[7] = algebraic[82]-algebraic[59] algebraic[29] = constants[39]-states[8] algebraic[83] = (constants[44]*algebraic[67]*algebraic[29])/(constants[45]+algebraic[29]) algebraic[31] = (constants[46]*states[8])/(constants[47]+states[8]) rates[8] = algebraic[83]-algebraic[31] algebraic[33] = constants[50]-states[9] algebraic[84] = (constants[49]*constants[53]*algebraic[77]*algebraic[33])/(constants[54]+constants[49]*algebraic[33]) algebraic[35] = (constants[49]*constants[57]*constants[52]*states[9])/(constants[58]+constants[49]*states[9]) rates[9] = algebraic[84]-algebraic[35] algebraic[37] = constants[50]-states[10] algebraic[85] = (constants[49]*constants[53]*algebraic[77]*algebraic[37])/(constants[54]+constants[49]*algebraic[37]) algebraic[39] = (constants[49]*constants[55]*constants[51]*states[10])/(constants[56]+constants[49]*states[10]) rates[10] = algebraic[85]-algebraic[39] algebraic[87] = constants[9]*algebraic[67]*states[0] algebraic[22] = constants[10]*states[2] rates[0] = ((algebraic[11]-algebraic[79])+algebraic[22])-algebraic[87] rates[2] = algebraic[87]-algebraic[22] algebraic[34] = (1.00000/(constants[128]*constants[116]))*log(constants[64]/states[11]) algebraic[40] = constants[67]*(power(states[14], 3.00000))*states[15]*states[16]*(states[17]-algebraic[34]) algebraic[92] = 1.00000*exp(constants[84]*states[17]*constants[128])*(power(states[11], 3.00000))*constants[66] algebraic[93] = 1.00000*exp((constants[84]-1.00000)*states[17]*constants[128])*(power(constants[64], 3.00000))*states[13] algebraic[94] = (constants[80]/((power(constants[81], 3.00000)+power(constants[64], 3.00000))*(constants[82]+constants[66])*(1.00000+constants[83]*exp((constants[84]-1.00000)*states[17]*constants[128]))))*(algebraic[92]-algebraic[93]) algebraic[95] = 1.00000/(1.00000+0.124500*exp(-0.100000*states[17]*constants[128])+0.0365000*constants[130]*exp(-states[17]*constants[128])) algebraic[96] = (((constants[85]*algebraic[95])/(1.00000+power(constants[86]/states[11], 1.50000)))*constants[65])/(constants[65]+constants[87]) algebraic[99] = constants[91]*(states[17]-algebraic[34]) algebraic[100] = algebraic[40]+algebraic[99]+3.00000*algebraic[94]+3.00000*algebraic[96] rates[11] = (constants[123]*algebraic[100]*constants[62])/(constants[59]*constants[116]*constants[115]) algebraic[46] = (constants[119]*constants[77]*states[17]*constants[115]*constants[128]*(states[12]*exp(states[17]*constants[128])-constants[65]))/(exp(states[17]*constants[128])-1.00000) algebraic[32] = states[10]/constants[50] algebraic[48] = 0.500000*((0.400000*algebraic[32])/constants[126]+0.600000) algebraic[44] = (constants[119]*constants[76]*4.00000*states[17]*constants[115]*constants[128]*(0.00100000*exp(2.00000*states[17]*constants[128])-0.341000*constants[66]))/(exp(2.00000*states[17]*constants[128])-1.00000) algebraic[51] = algebraic[44]*constants[78]*algebraic[48]*(power(states[18], 4.00000))*states[20]*states[21]*states[22] algebraic[53] = (algebraic[46]/(1.00000+algebraic[51]/constants[79]))*constants[78]*algebraic[48]*(power(states[18], 4.00000))*states[20]*states[21]*states[22] algebraic[36] = (1.00000/(constants[128]*constants[117]))*log(constants[65]/states[12]) algebraic[58] = constants[68]*states[23]*(0.886000*states[24]+0.114000*states[25])*(states[17]-algebraic[36]) algebraic[60] = constants[69]*states[26]*states[27]*(states[17]-algebraic[36]) algebraic[78] = 1.02000/(1.00000+exp(0.238500*((states[17]-algebraic[36])-59.2150))) algebraic[86] = (0.491240*exp(0.0803200*((states[17]+5.47600)-algebraic[36]))+exp(0.0617500*((states[17]-algebraic[36])-594.310)))/(1.00000+exp(-0.514300*((states[17]-algebraic[36])+4.75300))) algebraic[88] = algebraic[78]/(algebraic[78]+algebraic[86]) algebraic[89] = constants[70]*(power(constants[65]/5.40000, 1.0/2))*algebraic[88]*(states[17]-algebraic[36]) algebraic[90] = 1.00000/(1.00000+exp((7.48800-states[17])/5.98000)) algebraic[91] = constants[71]*algebraic[90]*(states[17]-algebraic[36]) algebraic[101] = ((algebraic[58]+algebraic[60]+algebraic[89]+algebraic[91])-2.00000*algebraic[96])+algebraic[53] rates[12] = (constants[123]*algebraic[101]*constants[62])/(constants[59]*constants[117]*constants[115]) algebraic[97] = (constants[88]*states[13])/(constants[89]+states[13]) algebraic[38] = (1.00000/(constants[128]*constants[118]))*log(constants[66]/states[13]) algebraic[98] = constants[90]*(states[17]-algebraic[38]) algebraic[102] = (algebraic[51]+algebraic[98]+algebraic[97])-2.00000*algebraic[94] algebraic[104] = custom_piecewise([greater(voi , 59.1000) & less(voi , 59.5000), constants[121] , True, constants[120]]) algebraic[106] = custom_piecewise([greater(sin(2.00000* pi*voi) , 0.999900), 10.0000 , True, 0.00000]) algebraic[108] = custom_piecewise([equal(constants[0] , 0.00000), 0.00000 , equal(constants[0] , 1.00000), algebraic[106] , True, (algebraic[104]-states[17])/constants[122]]) rates[17] = -1000.00*((algebraic[102]+algebraic[101]+algebraic[100])-algebraic[108]) rates[29] = custom_piecewise([greater(-1000.00*((algebraic[102]+algebraic[101]+algebraic[100])-algebraic[108]) , 30000.0), 1.00000-10000.0*states[29] , True, 1.00000]) algebraic[30] = algebraic[28]/constants[37] algebraic[111] = (constants[95]*(1.00000+2.00000*algebraic[30]))/(1.00000+2.00000*constants[124]) algebraic[112] = (constants[94]*(power(states[13], 2.00000)))/(power(algebraic[111], 2.00000)+power(states[13], 2.00000)) algebraic[113] = (constants[94]*states[30])/constants[96] algebraic[114] = (states[30]-states[28])/constants[105] rates[30] = (algebraic[112]-algebraic[113])-(algebraic[114]*constants[61])/constants[60] algebraic[103] = states[29]+0.00200000 algebraic[105] = 1.00000-exp(-algebraic[103]/constants[97]) algebraic[107] = exp(-algebraic[103]/constants[98]) algebraic[109] = constants[99]/(1.00000*(1.00000+exp((algebraic[102]+5.00000)/0.900000))) algebraic[110] = algebraic[109]*algebraic[105]*algebraic[107]*(states[28]-states[13]) algebraic[116] = 1.00000/(1.00000+(constants[103]*constants[104])/(power(constants[104]+states[28], 2.00000))) rates[28] = algebraic[116]*(algebraic[114]-algebraic[110]) algebraic[115] = (constants[106]*constants[109])/(power(constants[109]+states[13], 2.00000)) algebraic[117] = (constants[107]*constants[110])/(power(constants[110]+states[13], 2.00000)) algebraic[119] = (constants[108]*constants[111])/(power(constants[111]+states[13], 2.00000)) algebraic[120] = 1.00000/(1.00000+algebraic[117]+2.00000*algebraic[115]+algebraic[119]) rates[13] = algebraic[120]*(((constants[123]*algebraic[102]*constants[62])/(constants[59]*constants[118]*constants[115])-((algebraic[112]-algebraic[113])*constants[60])/constants[59])+(algebraic[110]*constants[61])/constants[59]) return(rates) def computeAlgebraic(constants, states, voi): algebraic = array([[0.0] * len(voi)] * sizeAlgebraic) states = array(states) voi = array(voi) algebraic[10] = 1.00000/(1.00000+exp((states[17]+87.5000)/10.3000)) algebraic[2] = (0.320000*(states[17]+47.1300))/(1.00000-exp(-0.100000*(states[17]+47.1300))) algebraic[12] = 0.0800000*exp(-states[17]/11.0000) algebraic[3] = custom_piecewise([greater_equal(states[17] , -40.0000), 0.00000 , True, 0.135000*exp((80.0000+states[17])/-6.80000)]) algebraic[13] = custom_piecewise([greater_equal(states[17] , -40.0000), 1.00000/(0.130000*(1.00000+exp(-(states[17]+10.6600)/11.1000))) , True, 3.56000*exp(0.0790000*states[17])+310000.*exp(0.350000*states[17])]) algebraic[4] = custom_piecewise([greater_equal(states[17] , -40.0000), 0.00000 , True, ((-127140.*exp(0.244400*states[17])-3.47400e-05*exp(-0.0439100*states[17]))*(states[17]+37.7800))/(1.00000+exp(0.311000*(states[17]+79.2300)))]) algebraic[14] = custom_piecewise([greater_equal(states[17] , -40.0000), (0.300000*exp(-2.57500e-07*states[17]))/(1.00000+exp(-0.100000*(states[17]+32.0000))) , True, (0.121200*exp(-0.0105200*states[17]))/(1.00000+exp(-0.137800*(states[17]+40.1400)))]) algebraic[5] = 400.000*exp((states[17]+2.00000)/10.0000) algebraic[15] = 50.0000*exp((-1.00000*(states[17]+2.00000))/13.0000) algebraic[1] = states[9]/constants[50] algebraic[17] = constants[72]*((0.375000*algebraic[1])/constants[125]+0.625000) algebraic[6] = 1.00000/(1.00000+exp((states[17]+55.0000)/7.50000))+0.100000/(1.00000+exp((-states[17]+21.0000)/6.00000)) algebraic[16] = 0.0200000+0.300000/(1.00000+exp((states[17]+30.0000)/9.50000)) algebraic[7] = 1.00000/(1.00000+exp((states[17]+10.6000)/-11.4200)) algebraic[18] = 1.00000/(45.1600*exp(0.0357700*(states[17]+50.0000))+98.9000*exp(-0.100000*(states[17]+38.0000))) algebraic[8] = 1.00000/(1.00000+exp((states[17]+43.5000)/6.88410)) algebraic[19] = 0.350000*exp(-1.00000*(power((states[17]+70.0000)/15.0000, 2.00000)))+0.0350000 algebraic[20] = 3.70000*exp(-1.00000*(power((states[17]+70.0000)/30.0000, 2.00000)))+0.0350000 algebraic[9] = 1.00000/(1.00000+exp((states[17]+11.5000)/-11.8200)) algebraic[21] = 10.0000/(45.1600*exp(0.0357700*(states[17]+50.0000))+98.9000*exp(-0.100000*(states[17]+38.0000))) algebraic[24] = constants[12]*states[3] algebraic[26] = constants[13]*states[4]*states[5] algebraic[23] = constants[74]*states[13] algebraic[25] = algebraic[23]*(power(1.00000-states[18], 4.00000)+2.00000*states[18]*(power(1.00000-states[18], 3.00000))+4.00000*(power(states[18], 2.00000))*(power(1.00000-states[18], 2.00000))+8.00000*(power(states[18], 3.00000))*(1.00000-states[18])+16.0000*(power(states[18], 4.00000))*(1.00000-algebraic[17]/constants[73])) algebraic[27] = constants[75]*(power(1.00000-states[19], 4.00000)+0.500000*states[19]*(power(1.00000-states[19], 3.00000))+0.250000*(power(states[19], 2.00000))*(power(1.00000-states[19], 2.00000))+0.125000*(power(states[19], 3.00000))*(1.00000-states[19])+0.0625000*(power(states[19], 4.00000))) algebraic[79] = constants[7]*(algebraic[61]+algebraic[62]) algebraic[11] = constants[8]*states[1] algebraic[80] = constants[11]*(algebraic[63]+algebraic[62]) algebraic[47] = (constants[20]*algebraic[43]*constants[15])/(constants[24]+constants[15]) algebraic[45] = (constants[19]*algebraic[42]*constants[15])/(constants[23]+constants[15]) algebraic[52] = (constants[21]*algebraic[50]*constants[15])/(constants[25]+constants[15]) algebraic[81] = (constants[22]*constants[131]*algebraic[68])/(constants[26]+algebraic[68]) algebraic[28] = constants[37]-states[7] algebraic[82] = (constants[40]*algebraic[67]*algebraic[28])/(constants[41]+algebraic[28]) algebraic[59] = (constants[42]*algebraic[56]*states[7])/(constants[43]+states[7]) algebraic[29] = constants[39]-states[8] algebraic[83] = (constants[44]*algebraic[67]*algebraic[29])/(constants[45]+algebraic[29]) algebraic[31] = (constants[46]*states[8])/(constants[47]+states[8]) algebraic[33] = constants[50]-states[9] algebraic[84] = (constants[49]*constants[53]*algebraic[77]*algebraic[33])/(constants[54]+constants[49]*algebraic[33]) algebraic[35] = (constants[49]*constants[57]*constants[52]*states[9])/(constants[58]+constants[49]*states[9]) algebraic[37] = constants[50]-states[10] algebraic[85] = (constants[49]*constants[53]*algebraic[77]*algebraic[37])/(constants[54]+constants[49]*algebraic[37]) algebraic[39] = (constants[49]*constants[55]*constants[51]*states[10])/(constants[56]+constants[49]*states[10]) algebraic[87] = constants[9]*algebraic[67]*states[0] algebraic[22] = constants[10]*states[2] algebraic[34] = (1.00000/(constants[128]*constants[116]))*log(constants[64]/states[11]) algebraic[40] = constants[67]*(power(states[14], 3.00000))*states[15]*states[16]*(states[17]-algebraic[34]) algebraic[92] = 1.00000*exp(constants[84]*states[17]*constants[128])*(power(states[11], 3.00000))*constants[66] algebraic[93] = 1.00000*exp((constants[84]-1.00000)*states[17]*constants[128])*(power(constants[64], 3.00000))*states[13] algebraic[94] = (constants[80]/((power(constants[81], 3.00000)+power(constants[64], 3.00000))*(constants[82]+constants[66])*(1.00000+constants[83]*exp((constants[84]-1.00000)*states[17]*constants[128]))))*(algebraic[92]-algebraic[93]) algebraic[95] = 1.00000/(1.00000+0.124500*exp(-0.100000*states[17]*constants[128])+0.0365000*constants[130]*exp(-states[17]*constants[128])) algebraic[96] = (((constants[85]*algebraic[95])/(1.00000+power(constants[86]/states[11], 1.50000)))*constants[65])/(constants[65]+constants[87]) algebraic[99] = constants[91]*(states[17]-algebraic[34]) algebraic[100] = algebraic[40]+algebraic[99]+3.00000*algebraic[94]+3.00000*algebraic[96] algebraic[46] = (constants[119]*constants[77]*states[17]*constants[115]*constants[128]*(states[12]*exp(states[17]*constants[128])-constants[65]))/(exp(states[17]*constants[128])-1.00000) algebraic[32] = states[10]/constants[50] algebraic[48] = 0.500000*((0.400000*algebraic[32])/constants[126]+0.600000) algebraic[44] = (constants[119]*constants[76]*4.00000*states[17]*constants[115]*constants[128]*(0.00100000*exp(2.00000*states[17]*constants[128])-0.341000*constants[66]))/(exp(2.00000*states[17]*constants[128])-1.00000) algebraic[51] = algebraic[44]*constants[78]*algebraic[48]*(power(states[18], 4.00000))*states[20]*states[21]*states[22] algebraic[53] = (algebraic[46]/(1.00000+algebraic[51]/constants[79]))*constants[78]*algebraic[48]*(power(states[18], 4.00000))*states[20]*states[21]*states[22] algebraic[36] = (1.00000/(constants[128]*constants[117]))*log(constants[65]/states[12]) algebraic[58] = constants[68]*states[23]*(0.886000*states[24]+0.114000*states[25])*(states[17]-algebraic[36]) algebraic[60] = constants[69]*states[26]*states[27]*(states[17]-algebraic[36]) algebraic[78] = 1.02000/(1.00000+exp(0.238500*((states[17]-algebraic[36])-59.2150))) algebraic[86] = (0.491240*exp(0.0803200*((states[17]+5.47600)-algebraic[36]))+exp(0.0617500*((states[17]-algebraic[36])-594.310)))/(1.00000+exp(-0.514300*((states[17]-algebraic[36])+4.75300))) algebraic[88] = algebraic[78]/(algebraic[78]+algebraic[86]) algebraic[89] = constants[70]*(power(constants[65]/5.40000, 1.0/2))*algebraic[88]*(states[17]-algebraic[36]) algebraic[90] = 1.00000/(1.00000+exp((7.48800-states[17])/5.98000)) algebraic[91] = constants[71]*algebraic[90]*(states[17]-algebraic[36]) algebraic[101] = ((algebraic[58]+algebraic[60]+algebraic[89]+algebraic[91])-2.00000*algebraic[96])+algebraic[53] algebraic[97] = (constants[88]*states[13])/(constants[89]+states[13]) algebraic[38] = (1.00000/(constants[128]*constants[118]))*log(constants[66]/states[13]) algebraic[98] = constants[90]*(states[17]-algebraic[38]) algebraic[102] = (algebraic[51]+algebraic[98]+algebraic[97])-2.00000*algebraic[94] algebraic[104] = custom_piecewise([greater(voi , 59.1000) & less(voi , 59.5000), constants[121] , True, constants[120]]) algebraic[106] = custom_piecewise([greater(sin(2.00000* pi*voi) , 0.999900), 10.0000 , True, 0.00000]) algebraic[108] = custom_piecewise([equal(constants[0] , 0.00000), 0.00000 , equal(constants[0] , 1.00000), algebraic[106] , True, (algebraic[104]-states[17])/constants[122]]) algebraic[30] = algebraic[28]/constants[37] algebraic[111] = (constants[95]*(1.00000+2.00000*algebraic[30]))/(1.00000+2.00000*constants[124]) algebraic[112] = (constants[94]*(power(states[13], 2.00000)))/(power(algebraic[111], 2.00000)+power(states[13], 2.00000)) algebraic[113] = (constants[94]*states[30])/constants[96] algebraic[114] = (states[30]-states[28])/constants[105] algebraic[103] = states[29]+0.00200000 algebraic[105] = 1.00000-exp(-algebraic[103]/constants[97]) algebraic[107] = exp(-algebraic[103]/constants[98]) algebraic[109] = constants[99]/(1.00000*(1.00000+exp((algebraic[102]+5.00000)/0.900000))) algebraic[110] = algebraic[109]*algebraic[105]*algebraic[107]*(states[28]-states[13]) algebraic[116] = 1.00000/(1.00000+(constants[103]*constants[104])/(power(constants[104]+states[28], 2.00000))) algebraic[115] = (constants[106]*constants[109])/(power(constants[109]+states[13], 2.00000)) algebraic[117] = (constants[107]*constants[110])/(power(constants[110]+states[13], 2.00000)) algebraic[119] = (constants[108]*constants[111])/(power(constants[111]+states[13], 2.00000)) algebraic[120] = 1.00000/(1.00000+algebraic[117]+2.00000*algebraic[115]+algebraic[119]) algebraic[0] = states[7]/constants[37] algebraic[57] = algebraic[51]+algebraic[53] algebraic[118] = 1000.00*(((states[28]+states[28]/algebraic[116])*constants[61])/constants[59]+(states[30]*constants[60])/constants[59]) return algebraic initialGuess0 = None def rootfind_0(voi, constants, rates, states, algebraic): """Calculate values of algebraic variables for DAE""" from scipy.optimize import fsolve global initialGuess0 if initialGuess0 is None: initialGuess0 = ones(3)*0.1 if not iterable(voi): soln = fsolve(residualSN_0, initialGuess0, args=(algebraic, voi, constants, rates, states), xtol=1E-6) initialGuess0 = soln constants[131] = soln[0] constants[132] = soln[1] constants[133] = soln[2] else: for (i,t) in enumerate(voi): soln = fsolve(residualSN_0, initialGuess0, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6) initialGuess0 = soln constants[131][i] = soln[0] constants[132][i] = soln[1] constants[133][i] = soln[2] def residualSN_0(algebraicCandidate, algebraic, voi, constants, rates, states): resid = array([0.0] * 3) constants[131] = algebraicCandidate[0] constants[132] = algebraicCandidate[1] constants[133] = algebraicCandidate[2] resid[0] = (constants[133]-(constants[131]*constants[132])/constants[29]) resid[1] = (constants[131]-(constants[16]-constants[133])) resid[2] = (constants[132]-(constants[17]-constants[133])) return resid initialGuess1 = None def rootfind_1(voi, constants, rates, states, algebraic): """Calculate values of algebraic variables for DAE""" from scipy.optimize import fsolve global initialGuess1 if initialGuess1 is None: initialGuess1 = ones(17)*0.1 if not iterable(voi): soln = fsolve(residualSN_1, initialGuess1, args=(algebraic, voi, constants, rates, states), xtol=1E-6) initialGuess1 = soln algebraic[61] = soln[0] algebraic[62] = soln[1] algebraic[63] = soln[2] algebraic[64] = soln[3] algebraic[65] = soln[4] algebraic[66] = soln[5] algebraic[67] = soln[6] algebraic[68] = soln[7] algebraic[69] = soln[8] algebraic[70] = soln[9] algebraic[71] = soln[10] algebraic[72] = soln[11] algebraic[73] = soln[12] algebraic[74] = soln[13] algebraic[75] = soln[14] algebraic[76] = soln[15] algebraic[77] = soln[16] else: for (i,t) in enumerate(voi): soln = fsolve(residualSN_1, initialGuess1, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6) initialGuess1 = soln algebraic[61][i] = soln[0] algebraic[62][i] = soln[1] algebraic[63][i] = soln[2] algebraic[64][i] = soln[3] algebraic[65][i] = soln[4] algebraic[66][i] = soln[5] algebraic[67][i] = soln[6] algebraic[68][i] = soln[7] algebraic[69][i] = soln[8] algebraic[70][i] = soln[9] algebraic[71][i] = soln[10] algebraic[72][i] = soln[11] algebraic[73][i] = soln[12] algebraic[74][i] = soln[13] algebraic[75][i] = soln[14] algebraic[76][i] = soln[15] algebraic[77][i] = soln[16] def residualSN_1(algebraicCandidate, algebraic, voi, constants, rates, states): resid = array([0.0] * 17) algebraic[61] = algebraicCandidate[0] algebraic[62] = algebraicCandidate[1] algebraic[63] = algebraicCandidate[2] algebraic[64] = algebraicCandidate[3] algebraic[65] = algebraicCandidate[4] algebraic[66] = algebraicCandidate[5] algebraic[67] = algebraicCandidate[6] algebraic[68] = algebraicCandidate[7] algebraic[69] = algebraicCandidate[8] algebraic[70] = algebraicCandidate[9] algebraic[71] = algebraicCandidate[10] algebraic[72] = algebraicCandidate[11] algebraic[73] = algebraicCandidate[12] algebraic[74] = algebraicCandidate[13] algebraic[75] = algebraicCandidate[14] algebraic[76] = algebraicCandidate[15] algebraic[77] = algebraicCandidate[16] resid[0] = (algebraic[61]-(algebraic[64]*algebraic[65])/constants[4]) resid[1] = (algebraic[62]-(algebraic[61]*algebraic[66])/constants[5]) resid[2] = (algebraic[63]-(algebraic[65]*algebraic[66])/constants[6]) resid[3] = (algebraic[64]-((constants[1]-algebraic[61])-algebraic[62])) resid[4] = (algebraic[65]-(((states[0]-algebraic[61])-algebraic[62])-algebraic[63])) resid[5] = (algebraic[66]-((constants[3]-algebraic[62])-algebraic[63])) resid[6] = (algebraic[69]-(constants[32]*constants[36])/(constants[36]+algebraic[67]+algebraic[77])) resid[7] = (algebraic[70]-(algebraic[67]/constants[35])*algebraic[67]*(1.00000+algebraic[69]/constants[36])) resid[8] = (algebraic[71]-algebraic[67]*(1.00000+algebraic[69]/constants[36])) resid[9] = (algebraic[72]-(algebraic[77]/constants[35])*algebraic[77]*(1.00000+algebraic[69]/constants[36])) resid[10] = (algebraic[73]-algebraic[77]*(1.00000+algebraic[69]/constants[36])) resid[11] = (algebraic[74]-(constants[33]/algebraic[68])*algebraic[70]) resid[12] = (algebraic[75]-(constants[33]/algebraic[68])*algebraic[72]) resid[13] = (algebraic[76]-((constants[33]*constants[34])/constants[35]+(constants[33]*algebraic[68])/constants[35]+(power(algebraic[68], 2.00000))/constants[35])) resid[14] = (algebraic[68]-((states[6]-(algebraic[74]+2.00000*algebraic[70]+2.00000*algebraic[71]))-(algebraic[75]+2.00000*algebraic[72]+2.00000*algebraic[73]))) resid[15] = (constants[112]-(2.00000*constants[30]*(power(algebraic[68], 2.00000))-algebraic[67]*(1.00000+algebraic[69]/constants[36])*(algebraic[76]*algebraic[67]+power(algebraic[68], 2.00000)))) resid[16] = (constants[113]-(2.00000*constants[31]*(power(algebraic[68], 2.00000))-algebraic[77]*(1.00000+algebraic[69]/constants[36])*(algebraic[76]*algebraic[77]+power(algebraic[68], 2.00000)))) return resid initialGuess2 = None def rootfind_2(voi, constants, rates, states, algebraic): """Calculate values of algebraic variables for DAE""" from scipy.optimize import fsolve global initialGuess2 if initialGuess2 is None: initialGuess2 = ones(3)*0.1 if not iterable(voi): soln = fsolve(residualSN_2, initialGuess2, args=(algebraic, voi, constants, rates, states), xtol=1E-6) initialGuess2 = soln algebraic[41] = soln[0] algebraic[42] = soln[1] algebraic[43] = soln[2] else: for (i,t) in enumerate(voi): soln = fsolve(residualSN_2, initialGuess2, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6) initialGuess2 = soln algebraic[41][i] = soln[0] algebraic[42][i] = soln[1] algebraic[43][i] = soln[2] def residualSN_2(algebraicCandidate, algebraic, voi, constants, rates, states): resid = array([0.0] * 3) algebraic[41] = algebraicCandidate[0] algebraic[42] = algebraicCandidate[1] algebraic[43] = algebraicCandidate[2] resid[0] = (algebraic[43]-(algebraic[41]*algebraic[42])/constants[27]) resid[1] = (algebraic[41]-(states[3]-algebraic[43])) resid[2] = (algebraic[42]-(constants[14]-algebraic[43])) return resid initialGuess3 = None def rootfind_3(voi, constants, rates, states, algebraic): """Calculate values of algebraic variables for DAE""" from scipy.optimize import fsolve global initialGuess3 if initialGuess3 is None: initialGuess3 = ones(2)*0.1 if not iterable(voi): soln = fsolve(residualSN_3, initialGuess3, args=(algebraic, voi, constants, rates, states), xtol=1E-6) initialGuess3 = soln algebraic[49] = soln[0] algebraic[50] = soln[1] else: for (i,t) in enumerate(voi): soln = fsolve(residualSN_3, initialGuess3, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6) initialGuess3 = soln algebraic[49][i] = soln[0] algebraic[50][i] = soln[1] def residualSN_3(algebraicCandidate, algebraic, voi, constants, rates, states): resid = array([0.0] * 2) algebraic[49] = algebraicCandidate[0] algebraic[50] = algebraicCandidate[1] resid[0] = (algebraic[50]-(algebraic[49]*algebraic[42])/constants[28]) resid[1] = (algebraic[49]-(constants[18]-algebraic[50])) return resid initialGuess4 = None def rootfind_4(voi, constants, rates, states, algebraic): """Calculate values of algebraic variables for DAE""" from scipy.optimize import fsolve global initialGuess4 if initialGuess4 is None: initialGuess4 = ones(3)*0.1 if not iterable(voi): soln = fsolve(residualSN_4, initialGuess4, args=(algebraic, voi, constants, rates, states), xtol=1E-6) initialGuess4 = soln algebraic[54] = soln[0] algebraic[55] = soln[1] algebraic[56] = soln[2] else: for (i,t) in enumerate(voi): soln = fsolve(residualSN_4, initialGuess4, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6) initialGuess4 = soln algebraic[54][i] = soln[0] algebraic[55][i] = soln[1] algebraic[56][i] = soln[2] def residualSN_4(algebraicCandidate, algebraic, voi, constants, rates, states): resid = array([0.0] * 3) algebraic[54] = algebraicCandidate[0] algebraic[55] = algebraicCandidate[1] algebraic[56] = algebraicCandidate[2] resid[0] = (algebraic[54]-(algebraic[55]*algebraic[56])/constants[48]) resid[1] = (algebraic[55]-(states[8]-algebraic[54])) resid[2] = (algebraic[56]-(constants[38]-algebraic[54])) return resid def custom_piecewise(cases): """Compute result of a piecewise function""" return select(cases[0::2],cases[1::2]) def solve_model(): """Solve model with ODE solver""" from scipy.integrate import ode # Initialise constants and state variables (init_states, constants) = initConsts() # Set timespan to solve over voi = linspace(0, 10, 500) # Construct ODE object to solve r = ode(computeRates) r.set_integrator('vode', method='bdf', atol=1e-06, rtol=1e-06, max_step=1) r.set_initial_value(init_states, voi[0]) r.set_f_params(constants) # Solve model states = array([[0.0] * len(voi)] * sizeStates) states[:,0] = init_states for (i,t) in enumerate(voi[1:]): if r.successful(): r.integrate(t) states[:,i+1] = r.y else: break # Compute algebraic variables algebraic = computeAlgebraic(constants, states, voi) return (voi, states, algebraic) def plot_model(voi, states, algebraic): """Plot variables against variable of integration""" import pylab (legend_states, legend_algebraic, legend_voi, legend_constants) = createLegends() pylab.figure(1) pylab.plot(voi,vstack((states,algebraic)).T) pylab.xlabel(legend_voi) pylab.legend(legend_states + legend_algebraic, loc='best') pylab.show() if __name__ == "__main__": (voi, states, algebraic) = solve_model() plot_model(voi, states, algebraic)