- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2021-12-06 08:25:06+13:00
- Desc:
- Adding units submodule manually
- Permanent Source URI:
- https://models.cellml.org/workspace/705/rawfile/9f1753b7ba6c67dc1d2964813a1ed5f3322a73ce/python/FCU_adenylylCyclase_old.py
# ----------------------------------------------------------------
# This python script is the export of FCU_adenylyl_cyclase.cellml made on 22 Oct 2021
# ----------------------------------------------------------------
# Size of variable arrays:
sizeAlgebraic = 141
sizeStates = 50
sizeConstants = 120
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_1a in component BG_parameters (fmol_per_sec)"
legend_constants[1] = "kappa_1b in component BG_parameters (fmol_per_sec)"
legend_constants[2] = "kappa_2a in component BG_parameters (fmol_per_sec)"
legend_constants[3] = "kappa_2b in component BG_parameters (fmol_per_sec)"
legend_constants[4] = "kappa_3a in component BG_parameters (fmol_per_sec)"
legend_constants[5] = "kappa_3b in component BG_parameters (fmol_per_sec)"
legend_constants[6] = "kappa_4a in component BG_parameters (fmol_per_sec)"
legend_constants[7] = "kappa_4b in component BG_parameters (fmol_per_sec)"
legend_constants[8] = "kappa_5 in component BG_parameters (fmol_per_sec)"
legend_constants[9] = "kappa_6 in component BG_parameters (fmol_per_sec)"
legend_constants[10] = "kappa_7 in component BG_parameters (fmol_per_sec)"
legend_constants[11] = "kappa_GiAC in component BG_parameters (fmol_per_sec)"
legend_constants[12] = "kappa_sig1_B1 in component BG_parameters (fmol_per_sec)"
legend_constants[13] = "kappa_sig2_B1 in component BG_parameters (fmol_per_sec)"
legend_constants[14] = "kappa_sig3_B1 in component BG_parameters (fmol_per_sec)"
legend_constants[15] = "kappa_sig4_B1 in component BG_parameters (fmol_per_sec)"
legend_constants[16] = "kappa_sig1_M2 in component BG_parameters (fmol_per_sec)"
legend_constants[17] = "kappa_sig2_M2 in component BG_parameters (fmol_per_sec)"
legend_constants[18] = "kappa_sig3_M2 in component BG_parameters (fmol_per_sec)"
legend_constants[19] = "kappa_sig4_M2 in component BG_parameters (fmol_per_sec)"
legend_constants[20] = "kappa_act1_Gs in component BG_parameters (fmol_per_sec)"
legend_constants[21] = "kappa_act2_Gs in component BG_parameters (fmol_per_sec)"
legend_constants[22] = "kappa_hyd_Gs in component BG_parameters (fmol_per_sec)"
legend_constants[23] = "kappa_reassoc_Gs in component BG_parameters (fmol_per_sec)"
legend_constants[24] = "kappa_PiPhos in component BG_parameters (fmol_per_sec)"
legend_constants[25] = "kappa_act1_Gi in component BG_parameters (fmol_per_sec)"
legend_constants[26] = "kappa_act2_Gi in component BG_parameters (fmol_per_sec)"
legend_constants[27] = "kappa_hyd_Gi in component BG_parameters (fmol_per_sec)"
legend_constants[28] = "kappa_reassoc_Gi in component BG_parameters (fmol_per_sec)"
legend_constants[29] = "K_ATP in component BG_parameters (per_fmol)"
legend_constants[30] = "K_cAMP in component BG_parameters (per_fmol)"
legend_constants[31] = "K_AC in component BG_parameters (per_fmol)"
legend_constants[32] = "K_AC_ATP in component BG_parameters (per_fmol)"
legend_constants[33] = "K_aGs_GTP_AC in component BG_parameters (per_fmol)"
legend_constants[34] = "K_aGs_GTP_AC_ATP in component BG_parameters (per_fmol)"
legend_constants[35] = "K_FSK_AC in component BG_parameters (per_fmol)"
legend_constants[36] = "K_FSK_AC_ATP in component BG_parameters (per_fmol)"
legend_constants[37] = "K_PDE in component BG_parameters (per_fmol)"
legend_constants[38] = "K_PDE_cAMP in component BG_parameters (per_fmol)"
legend_constants[39] = "K_five_AMP in component BG_parameters (per_fmol)"
legend_constants[40] = "K_IBMX in component BG_parameters (per_fmol)"
legend_constants[41] = "K_PDEinh in component BG_parameters (per_fmol)"
legend_constants[42] = "K_aGs_GTP in component BG_parameters (per_fmol)"
legend_constants[43] = "K_FSK in component BG_parameters (per_fmol)"
legend_constants[44] = "K_aGi_GTP in component BG_parameters (per_fmol)"
legend_constants[45] = "K_ACinh in component BG_parameters (per_fmol)"
legend_constants[46] = "K_PPi in component BG_parameters (per_fmol)"
legend_constants[47] = "K_L_B1 in component BG_parameters (per_fmol)"
legend_constants[48] = "K_R_B1 in component BG_parameters (per_fmol)"
legend_constants[49] = "K_Gs in component BG_parameters (per_fmol)"
legend_constants[50] = "K_LR_B1 in component BG_parameters (per_fmol)"
legend_constants[51] = "K_R_B1Gs in component BG_parameters (per_fmol)"
legend_constants[52] = "K_LR_B1Gs in component BG_parameters (per_fmol)"
legend_constants[53] = "K_L_M2 in component BG_parameters (per_fmol)"
legend_constants[54] = "K_R_M2 in component BG_parameters (per_fmol)"
legend_constants[55] = "K_Gi in component BG_parameters (per_fmol)"
legend_constants[56] = "K_LR_M2 in component BG_parameters (per_fmol)"
legend_constants[57] = "K_R_M2Gi in component BG_parameters (per_fmol)"
legend_constants[58] = "K_LR_M2Gi in component BG_parameters (per_fmol)"
legend_constants[59] = "K_beta_gamma_Gs in component BG_parameters (per_fmol)"
legend_constants[60] = "K_aGs_GDP in component BG_parameters (per_fmol)"
legend_constants[61] = "K_Pi in component BG_parameters (per_fmol)"
legend_constants[62] = "K_GTP in component BG_parameters (per_fmol)"
legend_constants[63] = "K_GDP in component BG_parameters (per_fmol)"
legend_constants[64] = "K_beta_gamma_Gi in component BG_parameters (per_fmol)"
legend_constants[65] = "K_aGi_GDP in component BG_parameters (per_fmol)"
legend_voi = "time in component environment (second)"
legend_constants[66] = "vol_myo in component environment (pL)"
legend_constants[67] = "freq in component environment (dimensionless)"
legend_constants[68] = "stimSt in component environment (second)"
legend_constants[69] = "stimSt2 in component environment (second)"
legend_constants[70] = "stimDur in component environment (second)"
legend_constants[71] = "tRamp in component environment (second)"
legend_constants[72] = "stimMag in component environment (fmol)"
legend_constants[73] = "stimHolding in component environment (fmol)"
legend_constants[119] = "m in component environment (fmol_per_sec)"
legend_constants[74] = "q_ATP_init in component environment (fmol)"
legend_constants[75] = "q_AC_init in component environment (fmol)"
legend_constants[76] = "q_cAMP_init in component environment (fmol)"
legend_constants[77] = "q_AC_ATP_init in component environment (fmol)"
legend_constants[78] = "q_FSK_init in component environment (fmol)"
legend_constants[79] = "q_FSK_AC_init in component environment (fmol)"
legend_constants[80] = "q_FSK_AC_ATP_init in component environment (fmol)"
legend_constants[81] = "q_aGs_GTP_init in component environment (fmol)"
legend_constants[82] = "q_aGs_GTP_AC_init in component environment (fmol)"
legend_constants[83] = "q_aGs_GTP_AC_ATP_init in component environment (fmol)"
legend_constants[84] = "q_PDE_init in component environment (fmol)"
legend_constants[85] = "q_PDEinh_init in component environment (fmol)"
legend_constants[86] = "q_PDE_cAMP_init in component environment (fmol)"
legend_constants[87] = "q_IBMX_init in component environment (fmol)"
legend_constants[88] = "q_five_AMP_init in component environment (fmol)"
legend_constants[89] = "q_aGi_GTP_init in component environment (fmol)"
legend_constants[90] = "q_ACinh_init in component environment (fmol)"
legend_constants[91] = "q_PPi_init in component environment (fmol)"
legend_constants[92] = "q_R_B1_init in component environment (fmol)"
legend_constants[93] = "q_Gs_init in component environment (fmol)"
legend_algebraic[0] = "q_L_B1_init in component environment (fmol)"
legend_constants[94] = "q_LR_B1_init in component environment (fmol)"
legend_constants[95] = "q_R_B1Gs_init in component environment (fmol)"
legend_constants[96] = "q_LR_B1Gs_init in component environment (fmol)"
legend_algebraic[1] = "q_L_M2_init in component environment (fmol)"
legend_constants[97] = "q_R_M2_init in component environment (fmol)"
legend_constants[98] = "q_Gi_init in component environment (fmol)"
legend_constants[99] = "q_LR_M2_init in component environment (fmol)"
legend_constants[100] = "q_R_M2Gi_init in component environment (fmol)"
legend_constants[101] = "q_LR_M2Gi_init in component environment (fmol)"
legend_constants[102] = "q_beta_gamma_Gs_init in component environment (fmol)"
legend_constants[103] = "q_aGs_GDP_init in component environment (fmol)"
legend_constants[104] = "q_Pi_init in component environment (fmol)"
legend_constants[105] = "q_beta_gamma_Gi_init in component environment (fmol)"
legend_constants[106] = "q_aGi_GDP_init in component environment (fmol)"
legend_constants[107] = "q_R_B1GsP_init in component environment (fmol)"
legend_constants[108] = "q_LR_B1GsP_init in component environment (fmol)"
legend_constants[109] = "q_GTP_init in component environment (fmol)"
legend_constants[110] = "q_GDP_init in component environment (fmol)"
legend_algebraic[35] = "L_B1_T in component environment (fmol)"
legend_algebraic[59] = "L_M2_T in component environment (fmol)"
legend_algebraic[36] = "R_B1_T in component environment (fmol)"
legend_algebraic[60] = "R_M2_T in component environment (fmol)"
legend_algebraic[73] = "Gs_T in component environment (fmol)"
legend_algebraic[94] = "Gi_T in component environment (fmol)"
legend_algebraic[16] = "adenosine_T in component environment (fmol)"
legend_algebraic[2] = "q_ATP in component environment (fmol)"
legend_algebraic[3] = "q_cAMP in component environment (fmol)"
legend_algebraic[4] = "q_AC in component environment (fmol)"
legend_algebraic[5] = "q_AC_ATP in component environment (fmol)"
legend_algebraic[6] = "q_aGs_GTP_AC in component environment (fmol)"
legend_algebraic[7] = "q_aGs_GTP_AC_ATP in component environment (fmol)"
legend_algebraic[8] = "q_FSK_AC in component environment (fmol)"
legend_algebraic[9] = "q_FSK_AC_ATP in component environment (fmol)"
legend_algebraic[10] = "q_PDE in component environment (fmol)"
legend_algebraic[13] = "q_PDE_cAMP in component environment (fmol)"
legend_algebraic[14] = "q_five_AMP in component environment (fmol)"
legend_algebraic[17] = "q_IBMX in component environment (fmol)"
legend_algebraic[18] = "q_PDEinh in component environment (fmol)"
legend_algebraic[19] = "q_aGs_GTP in component environment (fmol)"
legend_algebraic[20] = "q_FSK in component environment (fmol)"
legend_algebraic[22] = "q_aGi_GTP in component environment (fmol)"
legend_algebraic[23] = "q_ACinh in component environment (fmol)"
legend_algebraic[25] = "q_PPi in component environment (fmol)"
legend_algebraic[24] = "q_L_B1 in component environment (fmol)"
legend_algebraic[26] = "q_R_B1 in component environment (fmol)"
legend_algebraic[27] = "q_Gs in component environment (fmol)"
legend_algebraic[29] = "q_LR_B1 in component environment (fmol)"
legend_algebraic[31] = "q_R_B1Gs in component environment (fmol)"
legend_algebraic[33] = "q_LR_B1Gs in component environment (fmol)"
legend_algebraic[37] = "q_L_M2 in component environment (fmol)"
legend_algebraic[44] = "q_R_M2 in component environment (fmol)"
legend_algebraic[47] = "q_Gi in component environment (fmol)"
legend_algebraic[50] = "q_LR_M2 in component environment (fmol)"
legend_algebraic[53] = "q_R_M2Gi in component environment (fmol)"
legend_algebraic[56] = "q_LR_M2Gi in component environment (fmol)"
legend_algebraic[61] = "q_beta_gamma_Gs in component environment (fmol)"
legend_algebraic[69] = "q_aGs_GDP in component environment (fmol)"
legend_algebraic[74] = "q_Pi in component environment (fmol)"
legend_algebraic[79] = "q_GTP in component environment (fmol)"
legend_algebraic[83] = "q_GDP in component environment (fmol)"
legend_algebraic[86] = "q_beta_gamma_Gi in component environment (fmol)"
legend_algebraic[93] = "q_aGi_GDP in component environment (fmol)"
legend_states[0] = "q_ATP in component cAMP (fmol)"
legend_states[1] = "q_cAMP in component cAMP (fmol)"
legend_states[2] = "q_AC in component cAMP (fmol)"
legend_states[3] = "q_AC_ATP in component cAMP (fmol)"
legend_states[4] = "q_aGs_GTP_AC in component cAMP (fmol)"
legend_states[5] = "q_aGs_GTP_AC_ATP in component cAMP (fmol)"
legend_states[6] = "q_FSK_AC in component cAMP (fmol)"
legend_states[7] = "q_FSK_AC_ATP in component cAMP (fmol)"
legend_states[8] = "q_PDE in component cAMP (fmol)"
legend_states[9] = "q_PDE_cAMP in component cAMP (fmol)"
legend_states[10] = "q_five_AMP in component cAMP (fmol)"
legend_states[11] = "q_IBMX in component cAMP (fmol)"
legend_states[12] = "q_PDEinh in component cAMP (fmol)"
legend_states[13] = "q_aGs_GTP in component cAMP (fmol)"
legend_states[14] = "q_FSK in component cAMP (fmol)"
legend_states[15] = "q_aGi_GTP in component cAMP (fmol)"
legend_states[16] = "q_ACinh in component cAMP (fmol)"
legend_states[17] = "q_PPi in component cAMP (fmol)"
legend_states[18] = "q_L_B1 in component LRGbinding_B1AR (fmol)"
legend_states[19] = "q_R_B1 in component LRGbinding_B1AR (fmol)"
legend_states[20] = "q_Gs in component LRGbinding_B1AR (fmol)"
legend_states[21] = "q_LR_B1 in component LRGbinding_B1AR (fmol)"
legend_states[22] = "q_R_B1Gs in component LRGbinding_B1AR (fmol)"
legend_states[23] = "q_LR_B1Gs in component LRGbinding_B1AR (fmol)"
legend_states[24] = "q_L_M2 in component LRGbinding_M2 (fmol)"
legend_states[25] = "q_R_M2 in component LRGbinding_M2 (fmol)"
legend_states[26] = "q_Gi in component LRGbinding_M2 (fmol)"
legend_states[27] = "q_LR_M2 in component LRGbinding_M2 (fmol)"
legend_states[28] = "q_R_M2Gi in component LRGbinding_M2 (fmol)"
legend_states[29] = "q_LR_M2Gi in component LRGbinding_M2 (fmol)"
legend_states[30] = "q_R_B1 in component GsProtein (fmol)"
legend_states[31] = "q_Gs in component GsProtein (fmol)"
legend_states[32] = "q_R_B1Gs in component GsProtein (fmol)"
legend_states[33] = "q_LR_B1 in component GsProtein (fmol)"
legend_states[34] = "q_LR_B1Gs in component GsProtein (fmol)"
legend_states[35] = "q_aGs_GTP in component GsProtein (fmol)"
legend_states[36] = "q_beta_gamma_Gs in component GsProtein (fmol)"
legend_states[37] = "q_aGs_GDP in component GsProtein (fmol)"
legend_states[38] = "q_Pi in component GsProtein (fmol)"
legend_states[39] = "q_GTP in component GsProtein (fmol)"
legend_states[40] = "q_GDP in component GsProtein (fmol)"
legend_states[41] = "q_R_M2 in component GiProtein (fmol)"
legend_states[42] = "q_Gi in component GiProtein (fmol)"
legend_states[43] = "q_R_M2Gi in component GiProtein (fmol)"
legend_states[44] = "q_LR_M2 in component GiProtein (fmol)"
legend_states[45] = "q_LR_M2Gi in component GiProtein (fmol)"
legend_states[46] = "q_aGi_GTP in component GiProtein (fmol)"
legend_states[47] = "q_beta_gamma_Gi in component GiProtein (fmol)"
legend_states[48] = "q_aGi_GDP in component GiProtein (fmol)"
legend_states[49] = "q_Pi in component GiProtein (fmol)"
legend_constants[111] = "q_GTP in component GiProtein (fmol)"
legend_constants[112] = "q_GDP in component GiProtein (fmol)"
legend_constants[113] = "R in component constants (J_per_K_per_mol)"
legend_constants[114] = "T in component constants (kelvin)"
legend_constants[115] = "F in component constants (C_per_mol)"
legend_algebraic[100] = "v1a in component cAMP (fmol_per_sec)"
legend_algebraic[104] = "v1b in component cAMP (fmol_per_sec)"
legend_algebraic[108] = "v2a in component cAMP (fmol_per_sec)"
legend_algebraic[111] = "v2b in component cAMP (fmol_per_sec)"
legend_algebraic[114] = "v3a in component cAMP (fmol_per_sec)"
legend_algebraic[117] = "v3b in component cAMP (fmol_per_sec)"
legend_algebraic[120] = "v4a in component cAMP (fmol_per_sec)"
legend_algebraic[124] = "v4b in component cAMP (fmol_per_sec)"
legend_algebraic[128] = "v5 in component cAMP (fmol_per_sec)"
legend_algebraic[121] = "v6 in component cAMP (fmol_per_sec)"
legend_algebraic[125] = "v7 in component cAMP (fmol_per_sec)"
legend_algebraic[129] = "vGiAC in component cAMP (fmol_per_sec)"
legend_algebraic[28] = "mu_ATP in component cAMP (J_per_mol)"
legend_algebraic[32] = "mu_AC in component cAMP (J_per_mol)"
legend_algebraic[30] = "mu_cAMP in component cAMP (J_per_mol)"
legend_algebraic[34] = "mu_AC_ATP in component cAMP (J_per_mol)"
legend_algebraic[84] = "mu_FSK in component cAMP (J_per_mol)"
legend_algebraic[48] = "mu_FSK_AC in component cAMP (J_per_mol)"
legend_algebraic[51] = "mu_FSK_AC_ATP in component cAMP (J_per_mol)"
legend_algebraic[80] = "mu_aGs_GTP in component cAMP (J_per_mol)"
legend_algebraic[38] = "mu_aGs_GTP_AC in component cAMP (J_per_mol)"
legend_algebraic[45] = "mu_aGs_GTP_AC_ATP in component cAMP (J_per_mol)"
legend_algebraic[54] = "mu_PDE in component cAMP (J_per_mol)"
legend_algebraic[75] = "mu_PDEinh in component cAMP (J_per_mol)"
legend_algebraic[57] = "mu_PDE_cAMP in component cAMP (J_per_mol)"
legend_algebraic[70] = "mu_IBMX in component cAMP (J_per_mol)"
legend_algebraic[62] = "mu_five_AMP in component cAMP (J_per_mol)"
legend_algebraic[87] = "mu_aGi_GTP in component cAMP (J_per_mol)"
legend_algebraic[90] = "mu_ACinh in component cAMP (J_per_mol)"
legend_algebraic[95] = "mu_PPi in component cAMP (J_per_mol)"
legend_constants[116] = "vol in component cAMP (pL)"
legend_algebraic[11] = "ATP_T in component cAMP (fmol)"
legend_algebraic[12] = "AC_T in component cAMP (fmol)"
legend_algebraic[21] = "Gs_T in component cAMP (fmol)"
legend_algebraic[15] = "cAMP_T in component cAMP (fmol)"
legend_algebraic[39] = "mu_L_B1 in component LRGbinding_B1AR (J_per_mol)"
legend_algebraic[46] = "mu_R_B1 in component LRGbinding_B1AR (J_per_mol)"
legend_algebraic[49] = "mu_Gs in component LRGbinding_B1AR (J_per_mol)"
legend_algebraic[52] = "mu_LR_B1 in component LRGbinding_B1AR (J_per_mol)"
legend_algebraic[55] = "mu_R_B1Gs in component LRGbinding_B1AR (J_per_mol)"
legend_algebraic[58] = "mu_LR_B1Gs in component LRGbinding_B1AR (J_per_mol)"
legend_algebraic[63] = "vsig1_B1 in component LRGbinding_B1AR (fmol_per_sec)"
legend_algebraic[71] = "vsig2_B1 in component LRGbinding_B1AR (fmol_per_sec)"
legend_algebraic[76] = "vsig3_B1 in component LRGbinding_B1AR (fmol_per_sec)"
legend_algebraic[81] = "vsig4_B1 in component LRGbinding_B1AR (fmol_per_sec)"
legend_constants[117] = "vol in component LRGbinding_B1AR (pL)"
legend_algebraic[40] = "L_T in component LRGbinding_B1AR (fmol)"
legend_algebraic[41] = "R_T in component LRGbinding_B1AR (fmol)"
legend_algebraic[42] = "Gs_T in component LRGbinding_B1AR (fmol)"
legend_algebraic[64] = "mu_L_M2 in component LRGbinding_M2 (J_per_mol)"
legend_algebraic[72] = "mu_R_M2 in component LRGbinding_M2 (J_per_mol)"
legend_algebraic[77] = "mu_Gi in component LRGbinding_M2 (J_per_mol)"
legend_algebraic[82] = "mu_LR_M2 in component LRGbinding_M2 (J_per_mol)"
legend_algebraic[85] = "mu_R_M2Gi in component LRGbinding_M2 (J_per_mol)"
legend_algebraic[88] = "mu_LR_M2Gi in component LRGbinding_M2 (J_per_mol)"
legend_algebraic[91] = "vsig1_M2 in component LRGbinding_M2 (fmol_per_sec)"
legend_algebraic[96] = "vsig2_M2 in component LRGbinding_M2 (fmol_per_sec)"
legend_algebraic[101] = "vsig3_M2 in component LRGbinding_M2 (fmol_per_sec)"
legend_algebraic[105] = "vsig4_M2 in component LRGbinding_M2 (fmol_per_sec)"
legend_constants[118] = "vol in component LRGbinding_M2 (pL)"
legend_algebraic[65] = "L_T in component LRGbinding_M2 (fmol)"
legend_algebraic[66] = "R_T in component LRGbinding_M2 (fmol)"
legend_algebraic[67] = "Gi_T in component LRGbinding_M2 (fmol)"
legend_algebraic[130] = "vact1_Gs in component GsProtein (fmol_per_sec)"
legend_algebraic[132] = "vact2_Gs in component GsProtein (fmol_per_sec)"
legend_algebraic[134] = "vhyd_Gs in component GsProtein (fmol_per_sec)"
legend_algebraic[136] = "vreassoc_Gs in component GsProtein (fmol_per_sec)"
legend_algebraic[137] = "vPiPhos in component GsProtein (fmol_per_sec)"
legend_algebraic[89] = "mu_R_B1 in component GsProtein (J_per_mol)"
legend_algebraic[92] = "mu_Gs in component GsProtein (J_per_mol)"
legend_algebraic[97] = "mu_R_B1Gs in component GsProtein (J_per_mol)"
legend_algebraic[102] = "mu_LR_B1 in component GsProtein (J_per_mol)"
legend_algebraic[106] = "mu_LR_B1Gs in component GsProtein (J_per_mol)"
legend_algebraic[109] = "mu_aGs_GTP in component GsProtein (J_per_mol)"
legend_algebraic[112] = "mu_beta_gamma_Gs in component GsProtein (J_per_mol)"
legend_algebraic[115] = "mu_aGs_GDP in component GsProtein (J_per_mol)"
legend_algebraic[118] = "mu_Pi in component GsProtein (J_per_mol)"
legend_algebraic[122] = "mu_GTP in component GsProtein (J_per_mol)"
legend_algebraic[126] = "mu_GDP in component GsProtein (J_per_mol)"
legend_algebraic[43] = "R_T in component GsProtein (fmol)"
legend_algebraic[78] = "Gs_T in component GsProtein (fmol)"
legend_algebraic[135] = "vact1_Gi in component GiProtein (fmol_per_sec)"
legend_algebraic[138] = "vact2_Gi in component GiProtein (fmol_per_sec)"
legend_algebraic[139] = "vhyd_Gi in component GiProtein (fmol_per_sec)"
legend_algebraic[140] = "vreassoc_Gi in component GiProtein (fmol_per_sec)"
legend_algebraic[98] = "mu_R_M2 in component GiProtein (J_per_mol)"
legend_algebraic[103] = "mu_Gi in component GiProtein (J_per_mol)"
legend_algebraic[107] = "mu_R_M2Gi in component GiProtein (J_per_mol)"
legend_algebraic[110] = "mu_LR_M2 in component GiProtein (J_per_mol)"
legend_algebraic[113] = "mu_LR_M2Gi in component GiProtein (J_per_mol)"
legend_algebraic[116] = "mu_aGi_GTP in component GiProtein (J_per_mol)"
legend_algebraic[119] = "mu_beta_gamma_Gi in component GiProtein (J_per_mol)"
legend_algebraic[123] = "mu_aGi_GDP in component GiProtein (J_per_mol)"
legend_algebraic[127] = "mu_Pi in component GiProtein (J_per_mol)"
legend_algebraic[131] = "mu_GTP in component GiProtein (J_per_mol)"
legend_algebraic[133] = "mu_GDP in component GiProtein (J_per_mol)"
legend_algebraic[68] = "R_T in component GiProtein (fmol)"
legend_algebraic[99] = "Gi_T in component GiProtein (fmol)"
legend_rates[0] = "d/dt q_ATP in component cAMP (fmol)"
legend_rates[2] = "d/dt q_AC in component cAMP (fmol)"
legend_rates[3] = "d/dt q_AC_ATP in component cAMP (fmol)"
legend_rates[1] = "d/dt q_cAMP in component cAMP (fmol)"
legend_rates[14] = "d/dt q_FSK in component cAMP (fmol)"
legend_rates[6] = "d/dt q_FSK_AC in component cAMP (fmol)"
legend_rates[7] = "d/dt q_FSK_AC_ATP in component cAMP (fmol)"
legend_rates[13] = "d/dt q_aGs_GTP in component cAMP (fmol)"
legend_rates[4] = "d/dt q_aGs_GTP_AC in component cAMP (fmol)"
legend_rates[5] = "d/dt q_aGs_GTP_AC_ATP in component cAMP (fmol)"
legend_rates[9] = "d/dt q_PDE_cAMP in component cAMP (fmol)"
legend_rates[8] = "d/dt q_PDE in component cAMP (fmol)"
legend_rates[11] = "d/dt q_IBMX in component cAMP (fmol)"
legend_rates[12] = "d/dt q_PDEinh in component cAMP (fmol)"
legend_rates[10] = "d/dt q_five_AMP in component cAMP (fmol)"
legend_rates[15] = "d/dt q_aGi_GTP in component cAMP (fmol)"
legend_rates[16] = "d/dt q_ACinh in component cAMP (fmol)"
legend_rates[17] = "d/dt q_PPi in component cAMP (fmol)"
legend_rates[18] = "d/dt q_L_B1 in component LRGbinding_B1AR (fmol)"
legend_rates[19] = "d/dt q_R_B1 in component LRGbinding_B1AR (fmol)"
legend_rates[20] = "d/dt q_Gs in component LRGbinding_B1AR (fmol)"
legend_rates[21] = "d/dt q_LR_B1 in component LRGbinding_B1AR (fmol)"
legend_rates[22] = "d/dt q_R_B1Gs in component LRGbinding_B1AR (fmol)"
legend_rates[23] = "d/dt q_LR_B1Gs in component LRGbinding_B1AR (fmol)"
legend_rates[24] = "d/dt q_L_M2 in component LRGbinding_M2 (fmol)"
legend_rates[25] = "d/dt q_R_M2 in component LRGbinding_M2 (fmol)"
legend_rates[26] = "d/dt q_Gi in component LRGbinding_M2 (fmol)"
legend_rates[27] = "d/dt q_LR_M2 in component LRGbinding_M2 (fmol)"
legend_rates[28] = "d/dt q_R_M2Gi in component LRGbinding_M2 (fmol)"
legend_rates[29] = "d/dt q_LR_M2Gi in component LRGbinding_M2 (fmol)"
legend_rates[30] = "d/dt q_R_B1 in component GsProtein (fmol)"
legend_rates[32] = "d/dt q_R_B1Gs in component GsProtein (fmol)"
legend_rates[31] = "d/dt q_Gs in component GsProtein (fmol)"
legend_rates[33] = "d/dt q_LR_B1 in component GsProtein (fmol)"
legend_rates[34] = "d/dt q_LR_B1Gs in component GsProtein (fmol)"
legend_rates[35] = "d/dt q_aGs_GTP in component GsProtein (fmol)"
legend_rates[36] = "d/dt q_beta_gamma_Gs in component GsProtein (fmol)"
legend_rates[37] = "d/dt q_aGs_GDP in component GsProtein (fmol)"
legend_rates[38] = "d/dt q_Pi in component GsProtein (fmol)"
legend_rates[39] = "d/dt q_GTP in component GsProtein (fmol)"
legend_rates[40] = "d/dt q_GDP in component GsProtein (fmol)"
legend_rates[41] = "d/dt q_R_M2 in component GiProtein (fmol)"
legend_rates[43] = "d/dt q_R_M2Gi in component GiProtein (fmol)"
legend_rates[42] = "d/dt q_Gi in component GiProtein (fmol)"
legend_rates[44] = "d/dt q_LR_M2 in component GiProtein (fmol)"
legend_rates[45] = "d/dt q_LR_M2Gi in component GiProtein (fmol)"
legend_rates[46] = "d/dt q_aGi_GTP in component GiProtein (fmol)"
legend_rates[47] = "d/dt q_beta_gamma_Gi in component GiProtein (fmol)"
legend_rates[48] = "d/dt q_aGi_GDP in component GiProtein (fmol)"
legend_rates[49] = "d/dt q_Pi in component GiProtein (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] = 1.19049e+06
constants[1] = 0.00248563
constants[2] = 49840.4
constants[3] = 0.13449
constants[4] = 5.81948e+07
constants[5] = 6.76684e-17
constants[6] = 6013.05
constants[7] = 0.23128
constants[8] = 220.24
constants[9] = 1173.07
constants[10] = 136970
constants[11] = 351.702
constants[12] = 1949.85
constants[13] = 14.6972
constants[14] = 4.76612
constants[15] = 48932.3
constants[16] = 3468.13
constants[17] = 20.1802
constants[18] = 82.6965
constants[19] = 23570
constants[20] = 0.363553
constants[21] = 0.472994
constants[22] = 0.0560685
constants[23] = 1563.43
constants[24] = 428.713
constants[25] = 0.10939
constants[26] = 0.00358651
constants[27] = 0.413631
constants[28] = 3802.06
constants[29] = 6.842e-05
constants[30] = 0.0154937
constants[31] = 0.964843
constants[32] = 2.33903
constants[33] = 2.4781
constants[34] = 1.83726
constants[35] = 0.0212235
constants[36] = 0.0429592
constants[37] = 0.907054
constants[38] = 0.628454
constants[39] = 0.0154937
constants[40] = 0.0141004
constants[41] = 13.1991
constants[42] = 0.186657
constants[43] = 1.45328e-05
constants[44] = 0.024903
constants[45] = 8.26545
constants[46] = 0.000128372
constants[47] = 0.0783988
constants[48] = 0.0220282
constants[49] = 1.31984
constants[50] = 9.01184
constants[51] = 73.3398
constants[52] = 56.3705
constants[53] = 0.11083
constants[54] = 0.0323495
constants[55] = 0.501288
constants[56] = 1.35667
constants[57] = 37.7835
constants[58] = 23.0482
constants[59] = 0.0048057
constants[60] = 0.0612439
constants[61] = 0.546852
constants[62] = 0.000340182
constants[63] = 0.0718472
constants[64] = 0.0141226
constants[65] = 0.00843467
constants[66] = 34.4
constants[67] = 500
constants[68] = 3.5e-4
constants[69] = 3e-4
constants[70] = 0.25e-4
constants[71] = 1.8e-4
constants[72] = 1e1
constants[73] = 1e-5
constants[74] = 190
constants[75] = 1.889E-03
constants[76] = 1e-18
constants[77] = 1e-18
constants[78] = 3.800E-05
constants[79] = 1e-18
constants[80] = 1e-18
constants[81] = 1e-18
constants[82] = 1e-18
constants[83] = 1e-18
constants[84] = 1.482E-03
constants[85] = 1e-18
constants[86] = 1e-18
constants[87] = 3.80E-02
constants[88] = 1e-18
constants[89] = 1e-18
constants[90] = 1e-18
constants[91] = 1e-18
constants[92] = 0.0004579000
constants[93] = 0.1455400000
constants[94] = 1e-18
constants[95] = 1e-18
constants[96] = 1e-18
constants[97] = 0.00072
constants[98] = 0.00836
constants[99] = 1e-18
constants[100] = 1e-18
constants[101] = 1e-18
constants[102] = 1e-18
constants[103] = 1e-18
constants[104] = 570
constants[105] = 1e-18
constants[106] = 1e-18
constants[107] = 1e-18
constants[108] = 1e-18
constants[109] = 2.2
constants[110] = 1.1
states[0] = 1e-16
states[1] = 1e-16
states[2] = 1e-16
states[3] = 1e-16
states[4] = 1e-16
states[5] = 1e-16
states[6] = 1e-16
states[7] = 1e-16
states[8] = 1e-16
states[9] = 1e-16
states[10] = 1e-16
states[11] = 1e-16
states[12] = 1e-16
states[13] = 1e-16
states[14] = 1e-16
states[15] = 1e-16
states[16] = 1e-16
states[17] = 1e-16
states[18] = 1e-18
states[19] = 1e-18
states[20] = 1e-18
states[21] = 1e-18
states[22] = 1e-18
states[23] = 1e-18
states[24] = 1e-16
states[25] = 1e-16
states[26] = 1e-16
states[27] = 1e-16
states[28] = 1e-16
states[29] = 1e-16
states[30] = 1e-16
states[31] = 1e-16
states[32] = 1e-16
states[33] = 1e-16
states[34] = 1e-16
states[35] = 1e-16
states[36] = 1e-16
states[37] = 1e-16
states[38] = 1e-16
states[39] = 1e-16
states[40] = 1e-16
states[41] = 1e-16
states[42] = 1e-16
states[43] = 1e-16
states[44] = 1e-16
states[45] = 1e-16
states[46] = 1e-16
states[47] = 1e-16
states[48] = 1e-16
states[49] = 1e-16
constants[111] = 1e-16
constants[112] = 1e-16
constants[113] = 8.31
constants[114] = 310
constants[115] = 96485
constants[116] = 38.0
constants[117] = 34.4
constants[118] = 34.4
constants[119] = constants[72]/constants[71]
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
algebraic[26] = constants[92]+states[19]+states[30]
algebraic[46] = constants[113]*constants[114]*log(constants[48]*algebraic[26])
algebraic[27] = constants[93]+states[20]+states[31]
algebraic[49] = constants[113]*constants[114]*log(constants[49]*algebraic[27])
algebraic[31] = constants[95]+states[22]+states[32]
algebraic[55] = constants[113]*constants[114]*log(constants[51]*algebraic[31])
algebraic[63] = constants[12]*exp((algebraic[46]+algebraic[49])/(constants[113]*constants[114]))-exp(algebraic[55]/(constants[113]*constants[114]))
algebraic[0] = custom_piecewise([less(voi , constants[68]) & greater(voi , constants[68]-constants[71]), constants[73]+constants[119]*((voi-constants[68])+constants[71]) , greater_equal(voi , constants[68]) & less(voi , constants[68]+constants[70]), constants[72]+constants[73] , less(voi , constants[68]+constants[71]+constants[70]) & greater_equal(voi , constants[68]+constants[70]), constants[73]+-constants[119]*(((voi-constants[68])-constants[71])-constants[70]) , True, constants[73]])
algebraic[24] = algebraic[0]+states[18]
algebraic[39] = constants[113]*constants[114]*log(constants[47]*algebraic[24])
algebraic[33] = constants[96]+states[23]+states[34]
algebraic[58] = constants[113]*constants[114]*log(constants[52]*algebraic[33])
algebraic[71] = constants[13]*exp((algebraic[55]+algebraic[39])/(constants[113]*constants[114]))-exp(algebraic[58]/(constants[113]*constants[114]))
rates[22] = algebraic[63]-algebraic[71]
algebraic[29] = constants[94]+states[21]+states[33]
algebraic[52] = constants[113]*constants[114]*log(constants[50]*algebraic[29])
algebraic[76] = constants[14]*exp((algebraic[52]+algebraic[49])/(constants[113]*constants[114]))-exp(algebraic[58]/(constants[113]*constants[114]))
rates[20] = -algebraic[63]-algebraic[76]
rates[23] = algebraic[71]+algebraic[76]
algebraic[81] = constants[15]*exp((algebraic[46]+algebraic[39])/(constants[113]*constants[114]))-exp(algebraic[52]/(constants[113]*constants[114]))
rates[18] = -algebraic[71]-algebraic[81]
rates[19] = -algebraic[63]-algebraic[81]
rates[21] = -algebraic[76]+algebraic[81]
algebraic[44] = constants[97]+states[25]+states[41]
algebraic[72] = constants[113]*constants[114]*log(constants[54]*algebraic[44])
algebraic[47] = constants[98]+states[26]+states[42]
algebraic[77] = constants[113]*constants[114]*log(constants[55]*algebraic[47])
algebraic[53] = constants[100]+states[28]+states[43]
algebraic[85] = constants[113]*constants[114]*log(constants[57]*algebraic[53])
algebraic[91] = constants[16]*exp((algebraic[72]+algebraic[77])/(constants[113]*constants[114]))-exp(algebraic[85]/(constants[113]*constants[114]))
algebraic[1] = custom_piecewise([less(voi , constants[69]) & greater(voi , constants[69]-constants[71]), constants[73]+constants[119]*((voi-constants[69])+constants[71]) , greater_equal(voi , constants[69]) & less(voi , constants[69]+constants[70]), constants[72]+constants[73] , less(voi , constants[69]+constants[71]+constants[70]) & greater_equal(voi , constants[69]+constants[70]), constants[73]+-constants[119]*(((voi-constants[69])-constants[71])-constants[70]) , True, constants[73]])
algebraic[37] = algebraic[1]+states[24]
algebraic[64] = constants[113]*constants[114]*log(constants[53]*algebraic[37])
algebraic[56] = constants[101]+states[29]+states[45]
algebraic[88] = constants[113]*constants[114]*log(constants[58]*algebraic[56])
algebraic[96] = constants[17]*exp((algebraic[85]+algebraic[64])/(constants[113]*constants[114]))-exp(algebraic[88]/(constants[113]*constants[114]))
rates[28] = algebraic[91]-algebraic[96]
algebraic[50] = constants[99]+states[27]+states[44]
algebraic[82] = constants[113]*constants[114]*log(constants[56]*algebraic[50])
algebraic[101] = constants[18]*exp((algebraic[82]+algebraic[77])/(constants[113]*constants[114]))-exp(algebraic[88]/(constants[113]*constants[114]))
rates[26] = -algebraic[91]-algebraic[101]
rates[29] = algebraic[96]+algebraic[101]
algebraic[2] = constants[74]+states[0]
algebraic[28] = constants[113]*constants[114]*log(constants[29]*algebraic[2])
algebraic[4] = constants[75]+states[2]
algebraic[32] = constants[113]*constants[114]*log(constants[31]*algebraic[4])
algebraic[5] = constants[77]+states[3]
algebraic[34] = constants[113]*constants[114]*log(constants[32]*algebraic[5])
algebraic[100] = constants[0]*(exp((algebraic[32]+algebraic[28])/(constants[113]*constants[114]))-exp(algebraic[34]/(constants[113]*constants[114])))
algebraic[3] = constants[76]+states[1]
algebraic[30] = constants[113]*constants[114]*log(constants[30]*algebraic[3])
algebraic[25] = constants[91]+states[17]
algebraic[95] = constants[113]*constants[114]*log(constants[46]*algebraic[25])
algebraic[104] = constants[1]*(exp(algebraic[34]/(constants[113]*constants[114]))-exp((algebraic[32]+algebraic[30]+algebraic[95])/(constants[113]*constants[114])))
rates[3] = algebraic[100]-algebraic[104]
algebraic[105] = constants[19]*exp((algebraic[72]+algebraic[64])/(constants[113]*constants[114]))-exp(algebraic[82]/(constants[113]*constants[114]))
rates[24] = -algebraic[96]-algebraic[105]
rates[25] = -algebraic[91]-algebraic[105]
rates[27] = -algebraic[101]+algebraic[105]
algebraic[6] = constants[82]+states[4]
algebraic[38] = constants[113]*constants[114]*log(constants[33]*algebraic[6])
algebraic[7] = constants[83]+states[5]
algebraic[45] = constants[113]*constants[114]*log(constants[34]*algebraic[7])
algebraic[108] = constants[2]*(exp((algebraic[38]+algebraic[28])/(constants[113]*constants[114]))-exp(algebraic[45]/(constants[113]*constants[114])))
algebraic[111] = constants[3]*(exp(algebraic[45]/(constants[113]*constants[114]))-exp((algebraic[38]+algebraic[30]+algebraic[95])/(constants[113]*constants[114])))
rates[5] = algebraic[108]-algebraic[111]
algebraic[8] = constants[79]+states[6]
algebraic[48] = constants[113]*constants[114]*log(constants[35]*algebraic[8])
algebraic[9] = constants[80]+states[7]
algebraic[51] = constants[113]*constants[114]*log(constants[36]*algebraic[9])
algebraic[114] = constants[4]*(exp((algebraic[48]+algebraic[28])/(constants[113]*constants[114]))-exp(algebraic[51]/(constants[113]*constants[114])))
rates[0] = (-algebraic[100]-algebraic[114])-algebraic[108]
algebraic[117] = constants[5]*(exp(algebraic[51]/(constants[113]*constants[114]))-exp((algebraic[48]+algebraic[30]+algebraic[95])/(constants[113]*constants[114])))
rates[7] = algebraic[114]-algebraic[117]
algebraic[10] = constants[84]+states[8]
algebraic[54] = constants[113]*constants[114]*log(constants[37]*algebraic[10])
algebraic[13] = constants[86]+states[9]
algebraic[57] = constants[113]*constants[114]*log(constants[38]*algebraic[13])
algebraic[120] = constants[6]*(exp((algebraic[54]+algebraic[30])/(constants[113]*constants[114]))-exp(algebraic[57]/(constants[113]*constants[114])))
rates[1] = (algebraic[104]+algebraic[117]+algebraic[111])-algebraic[120]
algebraic[19] = constants[81]+states[13]+states[35]
algebraic[80] = constants[113]*constants[114]*log(constants[42]*algebraic[19])
algebraic[121] = constants[9]*(exp((algebraic[32]+algebraic[80])/(constants[113]*constants[114]))-exp(algebraic[38]/(constants[113]*constants[114])))
rates[13] = -algebraic[121]
rates[4] = (algebraic[121]-algebraic[108])+algebraic[111]
algebraic[20] = constants[78]+states[14]
algebraic[84] = constants[113]*constants[114]*log(constants[43]*algebraic[20])
algebraic[125] = constants[10]*(exp((algebraic[84]+algebraic[32])/(constants[113]*constants[114]))-exp(algebraic[48]/(constants[113]*constants[114])))
rates[14] = -algebraic[125]
rates[6] = (algebraic[125]+algebraic[117])-algebraic[114]
algebraic[14] = constants[88]+states[10]
algebraic[62] = constants[113]*constants[114]*log(constants[39]*algebraic[14])
algebraic[124] = constants[7]*(exp(algebraic[57]/(constants[113]*constants[114]))-exp((algebraic[54]+algebraic[62])/(constants[113]*constants[114])))
rates[9] = algebraic[120]-algebraic[124]
rates[10] = algebraic[124]
algebraic[22] = constants[89]+states[15]+states[46]
algebraic[87] = constants[113]*constants[114]*log(constants[44]*algebraic[22])
algebraic[23] = constants[90]+states[16]
algebraic[90] = constants[113]*constants[114]*log(constants[45]*algebraic[23])
algebraic[129] = constants[11]*(exp((algebraic[32]+algebraic[87])/(constants[113]*constants[114]))-exp(algebraic[90]/(constants[113]*constants[114])))
rates[2] = (((algebraic[104]-algebraic[100])-algebraic[121])-algebraic[125])-algebraic[129]
algebraic[18] = constants[85]+states[12]
algebraic[75] = constants[113]*constants[114]*log(constants[41]*algebraic[18])
algebraic[17] = constants[87]+states[11]
algebraic[70] = constants[113]*constants[114]*log(constants[40]*algebraic[17])
algebraic[128] = constants[8]*(exp((algebraic[54]+algebraic[70])/(constants[113]*constants[114]))-exp(algebraic[75]/(constants[113]*constants[114])))
rates[8] = (algebraic[124]-algebraic[120])-algebraic[128]
rates[11] = -algebraic[128]
rates[12] = algebraic[128]
rates[15] = -algebraic[129]
rates[16] = algebraic[129]
rates[17] = algebraic[129]
algebraic[89] = constants[113]*constants[114]*log(constants[48]*algebraic[26])
algebraic[97] = constants[113]*constants[114]*log(constants[51]*algebraic[31])
algebraic[109] = constants[113]*constants[114]*log(constants[42]*algebraic[19])
algebraic[61] = constants[102]+states[36]
algebraic[112] = constants[113]*constants[114]*log(constants[59]*algebraic[61])
algebraic[79] = constants[109]+states[39]+constants[111]
algebraic[122] = constants[113]*constants[114]*log(constants[62]*algebraic[79])
algebraic[83] = constants[110]+states[40]+constants[112]
algebraic[126] = constants[113]*constants[114]*log(constants[63]*algebraic[83])
algebraic[130] = constants[20]*(exp((algebraic[97]+algebraic[122])/(constants[113]*constants[114]))-exp((algebraic[109]+algebraic[112]+algebraic[89]+algebraic[126])/(constants[113]*constants[114])))
rates[30] = algebraic[130]
rates[32] = -algebraic[130]
algebraic[102] = constants[113]*constants[114]*log(constants[50]*algebraic[29])
algebraic[106] = constants[113]*constants[114]*log(constants[52]*algebraic[33])
algebraic[132] = constants[21]*(exp((algebraic[106]+algebraic[122])/(constants[113]*constants[114]))-exp((algebraic[109]+algebraic[112]+algebraic[102]+algebraic[126])/(constants[113]*constants[114])))
rates[33] = algebraic[132]
rates[34] = -algebraic[132]
algebraic[69] = constants[103]+states[37]
algebraic[115] = constants[113]*constants[114]*log(constants[60]*algebraic[69])
algebraic[74] = constants[104]+states[38]+states[49]
algebraic[118] = constants[113]*constants[114]*log(constants[61]*algebraic[74])
algebraic[134] = constants[22]*(exp(algebraic[109]/(constants[113]*constants[114]))-exp((algebraic[115]+algebraic[118])/(constants[113]*constants[114])))
rates[35] = (algebraic[130]+algebraic[132])-algebraic[134]
algebraic[98] = constants[113]*constants[114]*log(constants[54]*algebraic[44])
algebraic[107] = constants[113]*constants[114]*log(constants[57]*algebraic[53])
algebraic[116] = constants[113]*constants[114]*log(constants[44]*algebraic[22])
algebraic[86] = constants[105]+states[47]
algebraic[119] = constants[113]*constants[114]*log(constants[64]*algebraic[86])
algebraic[131] = constants[113]*constants[114]*log(constants[62]*algebraic[79])
algebraic[133] = constants[113]*constants[114]*log(constants[63]*algebraic[83])
algebraic[135] = constants[25]*(exp((algebraic[107]+algebraic[131])/(constants[113]*constants[114]))-exp((algebraic[116]+algebraic[119]+algebraic[98]+algebraic[133])/(constants[113]*constants[114])))
rates[41] = algebraic[135]
rates[43] = -algebraic[135]
algebraic[92] = constants[113]*constants[114]*log(constants[49]*algebraic[27])
algebraic[136] = constants[23]*(exp((algebraic[115]+algebraic[112])/(constants[113]*constants[114]))-exp(algebraic[92]/(constants[113]*constants[114])))
rates[31] = algebraic[136]
rates[36] = (algebraic[130]+algebraic[132])-algebraic[136]
rates[37] = algebraic[134]-algebraic[136]
algebraic[137] = constants[24]*(exp((algebraic[126]+algebraic[118])/(constants[113]*constants[114]))-exp(algebraic[122]/(constants[113]*constants[114])))
rates[38] = -algebraic[137]+algebraic[134]
rates[39] = algebraic[137]
rates[40] = -algebraic[137]
algebraic[110] = constants[113]*constants[114]*log(constants[56]*algebraic[50])
algebraic[113] = constants[113]*constants[114]*log(constants[58]*algebraic[56])
algebraic[138] = constants[26]*(exp((algebraic[113]+algebraic[131])/(constants[113]*constants[114]))-exp((algebraic[116]+algebraic[119]+algebraic[110]+algebraic[133])/(constants[113]*constants[114])))
rates[44] = algebraic[138]
rates[45] = -algebraic[138]
algebraic[93] = constants[106]+states[48]
algebraic[123] = constants[113]*constants[114]*log(constants[65]*algebraic[93])
algebraic[127] = constants[113]*constants[114]*log(constants[61]*algebraic[74])
algebraic[139] = constants[27]*(exp(algebraic[116]/(constants[113]*constants[114]))-exp((algebraic[123]+algebraic[127])/(constants[113]*constants[114])))
rates[46] = (algebraic[135]+algebraic[138])-algebraic[139]
rates[49] = algebraic[139]
algebraic[103] = constants[113]*constants[114]*log(constants[55]*algebraic[47])
algebraic[140] = constants[28]*(exp((algebraic[123]+algebraic[119])/(constants[113]*constants[114]))-exp(algebraic[103]/(constants[113]*constants[114])))
rates[42] = algebraic[140]
rates[47] = (algebraic[135]+algebraic[138])-algebraic[140]
rates[48] = algebraic[139]-algebraic[140]
return(rates)
def computeAlgebraic(constants, states, voi):
algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
states = array(states)
voi = array(voi)
if False: #len(voi)>1:
# HACK: make any neative values positive
print('******************')
print('HARDCODING NUMBERS FOR Q_L_B1 AND Q_GS TO BE NON NEGATIVE')
print('******************')
algebraic[24] = [max(a,finfo(float).eps) for a in algebraic[24]] # q_L_B1
algebraic[27] = [max(a,finfo(float).eps) for a in algebraic[27]] # q_Gs
algebraic[26] = constants[92]+states[19]+states[30]
algebraic[46] = constants[113]*constants[114]*log(constants[48]*algebraic[26])
algebraic[27] = constants[93]+states[20]+states[31]
algebraic[49] = constants[113]*constants[114]*log(constants[49]*algebraic[27])
algebraic[31] = constants[95]+states[22]+states[32]
algebraic[55] = constants[113]*constants[114]*log(constants[51]*algebraic[31])
algebraic[63] = constants[12]*exp((algebraic[46]+algebraic[49])/(constants[113]*constants[114]))-exp(algebraic[55]/(constants[113]*constants[114]))
algebraic[0] = custom_piecewise([less(voi , constants[68]) & greater(voi , constants[68]-constants[71]), constants[73]+constants[119]*((voi-constants[68])+constants[71]) , greater_equal(voi , constants[68]) & less(voi , constants[68]+constants[70]), constants[72]+constants[73] , less(voi , constants[68]+constants[71]+constants[70]) & greater_equal(voi , constants[68]+constants[70]), constants[73]+-constants[119]*(((voi-constants[68])-constants[71])-constants[70]) , True, constants[73]])
algebraic[24] = algebraic[0]+states[18]
algebraic[39] = constants[113]*constants[114]*log(constants[47]*algebraic[24])
algebraic[33] = constants[96]+states[23]+states[34]
algebraic[58] = constants[113]*constants[114]*log(constants[52]*algebraic[33])
algebraic[71] = constants[13]*exp((algebraic[55]+algebraic[39])/(constants[113]*constants[114]))-exp(algebraic[58]/(constants[113]*constants[114]))
algebraic[29] = constants[94]+states[21]+states[33]
algebraic[52] = constants[113]*constants[114]*log(constants[50]*algebraic[29])
algebraic[76] = constants[14]*exp((algebraic[52]+algebraic[49])/(constants[113]*constants[114]))-exp(algebraic[58]/(constants[113]*constants[114]))
algebraic[81] = constants[15]*exp((algebraic[46]+algebraic[39])/(constants[113]*constants[114]))-exp(algebraic[52]/(constants[113]*constants[114]))
algebraic[44] = constants[97]+states[25]+states[41]
algebraic[72] = constants[113]*constants[114]*log(constants[54]*algebraic[44])
algebraic[47] = constants[98]+states[26]+states[42]
algebraic[77] = constants[113]*constants[114]*log(constants[55]*algebraic[47])
algebraic[53] = constants[100]+states[28]+states[43]
algebraic[85] = constants[113]*constants[114]*log(constants[57]*algebraic[53])
algebraic[91] = constants[16]*exp((algebraic[72]+algebraic[77])/(constants[113]*constants[114]))-exp(algebraic[85]/(constants[113]*constants[114]))
algebraic[1] = custom_piecewise([less(voi , constants[69]) & greater(voi , constants[69]-constants[71]), constants[73]+constants[119]*((voi-constants[69])+constants[71]) , greater_equal(voi , constants[69]) & less(voi , constants[69]+constants[70]), constants[72]+constants[73] , less(voi , constants[69]+constants[71]+constants[70]) & greater_equal(voi , constants[69]+constants[70]), constants[73]+-constants[119]*(((voi-constants[69])-constants[71])-constants[70]) , True, constants[73]])
algebraic[37] = algebraic[1]+states[24]
algebraic[64] = constants[113]*constants[114]*log(constants[53]*algebraic[37])
algebraic[56] = constants[101]+states[29]+states[45]
algebraic[88] = constants[113]*constants[114]*log(constants[58]*algebraic[56])
algebraic[96] = constants[17]*exp((algebraic[85]+algebraic[64])/(constants[113]*constants[114]))-exp(algebraic[88]/(constants[113]*constants[114]))
algebraic[50] = constants[99]+states[27]+states[44]
algebraic[82] = constants[113]*constants[114]*log(constants[56]*algebraic[50])
algebraic[101] = constants[18]*exp((algebraic[82]+algebraic[77])/(constants[113]*constants[114]))-exp(algebraic[88]/(constants[113]*constants[114]))
algebraic[2] = constants[74]+states[0]
algebraic[28] = constants[113]*constants[114]*log(constants[29]*algebraic[2])
algebraic[4] = constants[75]+states[2]
algebraic[32] = constants[113]*constants[114]*log(constants[31]*algebraic[4])
algebraic[5] = constants[77]+states[3]
algebraic[34] = constants[113]*constants[114]*log(constants[32]*algebraic[5])
algebraic[100] = constants[0]*(exp((algebraic[32]+algebraic[28])/(constants[113]*constants[114]))-exp(algebraic[34]/(constants[113]*constants[114])))
algebraic[3] = constants[76]+states[1]
algebraic[30] = constants[113]*constants[114]*log(constants[30]*algebraic[3])
algebraic[25] = constants[91]+states[17]
algebraic[95] = constants[113]*constants[114]*log(constants[46]*algebraic[25])
algebraic[104] = constants[1]*(exp(algebraic[34]/(constants[113]*constants[114]))-exp((algebraic[32]+algebraic[30]+algebraic[95])/(constants[113]*constants[114])))
algebraic[105] = constants[19]*exp((algebraic[72]+algebraic[64])/(constants[113]*constants[114]))-exp(algebraic[82]/(constants[113]*constants[114]))
algebraic[6] = constants[82]+states[4]
algebraic[38] = constants[113]*constants[114]*log(constants[33]*algebraic[6])
algebraic[7] = constants[83]+states[5]
algebraic[45] = constants[113]*constants[114]*log(constants[34]*algebraic[7])
algebraic[108] = constants[2]*(exp((algebraic[38]+algebraic[28])/(constants[113]*constants[114]))-exp(algebraic[45]/(constants[113]*constants[114])))
algebraic[111] = constants[3]*(exp(algebraic[45]/(constants[113]*constants[114]))-exp((algebraic[38]+algebraic[30]+algebraic[95])/(constants[113]*constants[114])))
algebraic[8] = constants[79]+states[6]
algebraic[48] = constants[113]*constants[114]*log(constants[35]*algebraic[8])
algebraic[9] = constants[80]+states[7]
algebraic[51] = constants[113]*constants[114]*log(constants[36]*algebraic[9])
algebraic[114] = constants[4]*(exp((algebraic[48]+algebraic[28])/(constants[113]*constants[114]))-exp(algebraic[51]/(constants[113]*constants[114])))
algebraic[117] = constants[5]*(exp(algebraic[51]/(constants[113]*constants[114]))-exp((algebraic[48]+algebraic[30]+algebraic[95])/(constants[113]*constants[114])))
algebraic[10] = constants[84]+states[8]
algebraic[54] = constants[113]*constants[114]*log(constants[37]*algebraic[10])
algebraic[13] = constants[86]+states[9]
algebraic[57] = constants[113]*constants[114]*log(constants[38]*algebraic[13])
algebraic[120] = constants[6]*(exp((algebraic[54]+algebraic[30])/(constants[113]*constants[114]))-exp(algebraic[57]/(constants[113]*constants[114])))
algebraic[19] = constants[81]+states[13]+states[35]
algebraic[80] = constants[113]*constants[114]*log(constants[42]*algebraic[19])
algebraic[121] = constants[9]*(exp((algebraic[32]+algebraic[80])/(constants[113]*constants[114]))-exp(algebraic[38]/(constants[113]*constants[114])))
algebraic[20] = constants[78]+states[14]
algebraic[84] = constants[113]*constants[114]*log(constants[43]*algebraic[20])
algebraic[125] = constants[10]*(exp((algebraic[84]+algebraic[32])/(constants[113]*constants[114]))-exp(algebraic[48]/(constants[113]*constants[114])))
algebraic[14] = constants[88]+states[10]
algebraic[62] = constants[113]*constants[114]*log(constants[39]*algebraic[14])
algebraic[124] = constants[7]*(exp(algebraic[57]/(constants[113]*constants[114]))-exp((algebraic[54]+algebraic[62])/(constants[113]*constants[114])))
algebraic[22] = constants[89]+states[15]+states[46]
algebraic[87] = constants[113]*constants[114]*log(constants[44]*algebraic[22])
algebraic[23] = constants[90]+states[16]
algebraic[90] = constants[113]*constants[114]*log(constants[45]*algebraic[23])
algebraic[129] = constants[11]*(exp((algebraic[32]+algebraic[87])/(constants[113]*constants[114]))-exp(algebraic[90]/(constants[113]*constants[114])))
algebraic[18] = constants[85]+states[12]
algebraic[75] = constants[113]*constants[114]*log(constants[41]*algebraic[18])
algebraic[17] = constants[87]+states[11]
algebraic[70] = constants[113]*constants[114]*log(constants[40]*algebraic[17])
algebraic[128] = constants[8]*(exp((algebraic[54]+algebraic[70])/(constants[113]*constants[114]))-exp(algebraic[75]/(constants[113]*constants[114])))
algebraic[89] = constants[113]*constants[114]*log(constants[48]*algebraic[26])
algebraic[97] = constants[113]*constants[114]*log(constants[51]*algebraic[31])
algebraic[109] = constants[113]*constants[114]*log(constants[42]*algebraic[19])
algebraic[61] = constants[102]+states[36]
algebraic[112] = constants[113]*constants[114]*log(constants[59]*algebraic[61])
algebraic[79] = constants[109]+states[39]+constants[111]
algebraic[122] = constants[113]*constants[114]*log(constants[62]*algebraic[79])
algebraic[83] = constants[110]+states[40]+constants[112]
algebraic[126] = constants[113]*constants[114]*log(constants[63]*algebraic[83])
algebraic[130] = constants[20]*(exp((algebraic[97]+algebraic[122])/(constants[113]*constants[114]))-exp((algebraic[109]+algebraic[112]+algebraic[89]+algebraic[126])/(constants[113]*constants[114])))
algebraic[102] = constants[113]*constants[114]*log(constants[50]*algebraic[29])
algebraic[106] = constants[113]*constants[114]*log(constants[52]*algebraic[33])
algebraic[132] = constants[21]*(exp((algebraic[106]+algebraic[122])/(constants[113]*constants[114]))-exp((algebraic[109]+algebraic[112]+algebraic[102]+algebraic[126])/(constants[113]*constants[114])))
algebraic[69] = constants[103]+states[37]
algebraic[115] = constants[113]*constants[114]*log(constants[60]*algebraic[69])
algebraic[74] = constants[104]+states[38]+states[49]
algebraic[118] = constants[113]*constants[114]*log(constants[61]*algebraic[74])
algebraic[134] = constants[22]*(exp(algebraic[109]/(constants[113]*constants[114]))-exp((algebraic[115]+algebraic[118])/(constants[113]*constants[114])))
algebraic[98] = constants[113]*constants[114]*log(constants[54]*algebraic[44])
algebraic[107] = constants[113]*constants[114]*log(constants[57]*algebraic[53])
algebraic[116] = constants[113]*constants[114]*log(constants[44]*algebraic[22])
algebraic[86] = constants[105]+states[47]
algebraic[119] = constants[113]*constants[114]*log(constants[64]*algebraic[86])
algebraic[131] = constants[113]*constants[114]*log(constants[62]*algebraic[79])
algebraic[133] = constants[113]*constants[114]*log(constants[63]*algebraic[83])
algebraic[135] = constants[25]*(exp((algebraic[107]+algebraic[131])/(constants[113]*constants[114]))-exp((algebraic[116]+algebraic[119]+algebraic[98]+algebraic[133])/(constants[113]*constants[114])))
algebraic[92] = constants[113]*constants[114]*log(constants[49]*algebraic[27])
algebraic[136] = constants[23]*(exp((algebraic[115]+algebraic[112])/(constants[113]*constants[114]))-exp(algebraic[92]/(constants[113]*constants[114])))
algebraic[137] = constants[24]*(exp((algebraic[126]+algebraic[118])/(constants[113]*constants[114]))-exp(algebraic[122]/(constants[113]*constants[114])))
algebraic[110] = constants[113]*constants[114]*log(constants[56]*algebraic[50])
algebraic[113] = constants[113]*constants[114]*log(constants[58]*algebraic[56])
algebraic[138] = constants[26]*(exp((algebraic[113]+algebraic[131])/(constants[113]*constants[114]))-exp((algebraic[116]+algebraic[119]+algebraic[110]+algebraic[133])/(constants[113]*constants[114])))
algebraic[93] = constants[106]+states[48]
algebraic[123] = constants[113]*constants[114]*log(constants[65]*algebraic[93])
algebraic[127] = constants[113]*constants[114]*log(constants[61]*algebraic[74])
algebraic[139] = constants[27]*(exp(algebraic[116]/(constants[113]*constants[114]))-exp((algebraic[123]+algebraic[127])/(constants[113]*constants[114])))
algebraic[103] = constants[113]*constants[114]*log(constants[55]*algebraic[47])
algebraic[140] = constants[28]*(exp((algebraic[123]+algebraic[119])/(constants[113]*constants[114]))-exp(algebraic[103]/(constants[113]*constants[114])))
algebraic[11] = algebraic[2]+algebraic[5]+algebraic[9]+algebraic[7]
algebraic[12] = algebraic[4]+algebraic[5]+algebraic[8]+algebraic[9]+algebraic[6]+algebraic[7]
algebraic[15] = algebraic[3]+algebraic[13]+states[10]
algebraic[16] = algebraic[3]+algebraic[13]+algebraic[14]+algebraic[2]+algebraic[5]+algebraic[7]+algebraic[9]
algebraic[21] = algebraic[19]+algebraic[6]+algebraic[7]
algebraic[35] = algebraic[24]+algebraic[33]+algebraic[29]
algebraic[36] = algebraic[26]+algebraic[31]+algebraic[29]+algebraic[33]
algebraic[40] = algebraic[24]+algebraic[29]+algebraic[33]
algebraic[41] = algebraic[26]+algebraic[29]+algebraic[31]+algebraic[33]
algebraic[42] = algebraic[27]+algebraic[31]+algebraic[33]
algebraic[43] = algebraic[26]+algebraic[31]+algebraic[29]+algebraic[33]
algebraic[59] = algebraic[37]+algebraic[56]+algebraic[50]
algebraic[60] = algebraic[44]+algebraic[53]+algebraic[50]+algebraic[56]
algebraic[65] = algebraic[37]+algebraic[50]+algebraic[56]
algebraic[66] = algebraic[44]+algebraic[50]+algebraic[53]+algebraic[56]
algebraic[67] = algebraic[47]+algebraic[53]+algebraic[56]
algebraic[68] = algebraic[44]+algebraic[53]+algebraic[50]+algebraic[56]
algebraic[73] = algebraic[27]+algebraic[31]+algebraic[33]+algebraic[19]+algebraic[69]
algebraic[78] = algebraic[27]+algebraic[31]+algebraic[33]+algebraic[19]+algebraic[69]
algebraic[94] = algebraic[47]+algebraic[53]+algebraic[56]+algebraic[22]+algebraic[93]
algebraic[99] = algebraic[47]+algebraic[53]+algebraic[56]+algebraic[22]+algebraic[93]
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_AC_init in component environment (fmol)',
# 'q_cAMP_init in component environment (fmol)',
# 'q_R_B1_init in component environment (fmol)',
'q_Gs_init in component environment (fmol)',
# 'q_R_M2_init in component environment (fmol)',
# 'q_Gi_init in component environment (fmol)'
]
[igs,_,_] = find_indices(glabels, legend_constants, legend_states, legend_algebraic);
Gmult = 5e1; #1e3;
(init_states, constants) = initConsts(Gmult,igs)
# Set timespan to solve over
voi = linspace(0, 12e-5, 150000)
# options for stimulus protocol
stimLabels = ["stimSt in component environment (second)",
"stimSt2 in component environment (second)",
"stimDur in component environment (second)",
"tRamp in component environment (second)",
"stimMag in component environment (fmol)",
"stimHolding in component environment (fmol)",
"freq in component environment (dimensionless)",
"m in component environment (fmol_per_sec)"]
sd = find_indices_from_labels(stimLabels, legend_constants)
constants[sd['freq']] = 1 #[0.25,0.25]; # per second # UNUSED
constants[sd['stimSt']] = 5e-5
constants[sd['stimSt2']] = 4e-5
constants[sd['stimHolding']] = 1e-9
constants[sd['stimMag']] = 1.888e1 #[1e6,1e6];
constants[sd['tRamp']] = 5e-5 # seconds for rAMP start
constants[sd['stimDur']] = 1e-5
if False:
constants[sd['stimMag']] = 0
# recalculate m
constants[sd['m']] = constants[sd['stimMag']]/constants[sd['tRamp']]
# Construct ODE object to solve
r = ode(computeRates)
r.set_integrator('vode', method='bdf', atol=1e-07, rtol=1e-07, max_step=1e-9)
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 = ['q_L_B1_init in component environment (fmol)',
'q_L_M2_init in component environment (fmol)',
'q_LR_B1Gs in component environment (fmol)',
'q_LR_M2Gi in component environment (fmol)',
'q_aGs_GTP in component environment (fmol)',
'q_aGi_GTP in component environment (fmol)',
'q_aGs_GTP_AC in component cAMP (fmol)',
'q_AC in component environment (fmol)',
'q_ACinh in component environment (fmol)',
'q_cAMP in component environment (fmol)',
'q_ATP in component environment (fmol)']
if True:
[_, i_st, i_alg] = find_indices(labels, legend_constants, legend_states, legend_algebraic)
plot_selected(i_st,voi,states,legend_voi,legend_states,('LR_B1\nGmult=%g' %(Gmult)),int(ceil(sqrt(len(i_st)))))
plot_selected(i_alg,voi,algebraic,legend_voi,legend_algebraic,('LR_B1\nGmult=%g' %(Gmult)),int(ceil(sqrt(len(i_alg)))))
# plot chemostat concentrations: they should be constant, or machine
# error order
labels = [
'L_B1_T in component environment (fmol)',
'L_M2_T in component environment (fmol)',
'R_B1_T in component environment (fmol)',
'Gs_T in component environment (fmol)',
'R_M2_T in component environment (fmol)',
'Gi_T in component environment (fmol)',
'adenosine_T in component environment (fmol)'
]
labels_init = [
'q_L_B1_init in component environment (fmol)',
'q_L_M2_init in component environment (fmol)',
'q_R_B1_init in component environment (fmol)',
'q_Gs_init in component environment (fmol)',
'q_R_M2_init in component environment (fmol)',
'q_Gi_init in component environment (fmol)',
'q_ATP_init in component environment (fmol)'
]
[_, _, i_alg] = find_indices(labels, legend_constants, legend_states, legend_algebraic)
[i_con, _, i_algstim] = find_indices(labels, legend_constants, legend_states, legend_algebraic)
plot_2per(i_algstim,i_con, i_alg,voi,constants,algebraic,legend_algebraic,legend_constants,('chemostats\nGmult=%g' %(Gmult)),3)
# debug R_B1_t
if False:
labels = ['q_R_B1 in component LRGbinding_B1AR (fmol)','q_R_B1 in component GsProtein (fmol)']
[_, i_st, _] = find_indices(labels, legend_constants, legend_states, legend_algebraic)
plot_selected(i_st,VOI,states,legend_VOI,LEGEND_STATES,('R_B1\nGmult=%g' %(Gmult)),int(ceil(sqrt(len(i_st)))))
# debug aGs_GTP and aGiGTP
if False:
labels = [
'q_R_B1GsP in component environment (fmol)',
'q_R_M2GiP in component environment (fmol)',
'q_aGs_GTP in component environment (fmol)',
'q_aGi_GTP in component environment (fmol)',
'vphos_R_B1Gs in component GsProtein (fmol_per_sec)',
'vphos_LR_B1Gs in component GsProtein (fmol_per_sec)',
'vphos_R_M2Gi in component GiProtein (fmol_per_sec)',
'vphos_LR_M2Gi in component GiProtein (fmol_per_sec)',
'vact1_Gs in component GsProtein (fmol_per_sec)',
'vact2_Gs in component GsProtein (fmol_per_sec)',
'vact1_Gi in component GiProtein (fmol_per_sec)',
'vact2_Gi in component GiProtein (fmol_per_sec)',
'vhyd_Gs in component GsProtein (fmol_per_sec)',
'vhyd_Gi in component GiProtein (fmol_per_sec)',
'q_aGs_GTP in component GsProtein (fmol)',
'q_aGi_GTP in component GiProtein (fmol)',
]
[_, i_st, i_alg] = find_indices(labels, legend_constants, legend_states, legend_algebraic)
plot_selected(i_alg,voi,algebraic,legend_voi,legend_algebraic,('ag_gtp\ngmult=%g' %(Gmult)),int(ceil(sqrt(len(i_alg)))))
plot_selected(i_st,voi,states,legend_voi,legend_states,('ag_GTP\nGmult=%g' %(Gmult)),int(ceil(sqrt(len(i_st)))))
# 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])
if False:
plt.figure(1)
plt.plot(voi,vstack((states,algebraic)).T)
plt.xlabel(legend_voi)
plt.legend(legend_states + legend_algebraic, loc='best')
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)
def plot_2per(i_alg0,id0,ids,x,y0,y,legend_y,legend_y_con,titlestr,ns):
istart = 30;
plt.figure()
if i_alg0:
idoffset = len(i_alg0)
for i in range(idoffset):
plt.subplot(ns,ns,i+1)
plt.plot(x[istart:-1], y[i_alg0[i],istart:-1],'x-')
plt.plot(x[istart:-1], y[ids[i],istart:-1])
plt.xlabel('time (s)')
diff = mean(abs(y[i_alg0[i],istart:-1] - y[ids[i],istart:-1]))
plt.title(legend_y[i_alg0[i]].split(' in')[0] +' avg abs error = ' + str(diff))
str1 = legend_y[i_alg0[i]].split(' ')[0]
str2 = legend_y[ids[i]].split(' ')[0]
plt.legend([str1,str2])
else:
idoffset = 0
for i, id in enumerate(ids[idoffset+1:-1]):
plt.subplot(ns,ns,i+1)
plt.plot(x[istart:-1], (len(x)-istart+1)*[y0[id0[i-idoffset]]],'x-')
plt.plot(x[istart:-1], y[ids,istart:-1])
diff = abs(mean(y[ids, istart:-1]) - y0[id0[i-idoffset]])
plt.xlabel('time (s)')
str1 = legend_y[ids].split(' ')[0]
plt.title(str1 + ' avg abs error = ' + str(diff))
plt.legend([legend_y_con[id0[i-idoffset]],legend_y[ids]])
# plt.tight_layout()
plt.subplots_adjust(wspace=0.2, hspace=0.4)
plt.suptitle(titlestr +' eps = ' +str(finfo(float).eps))
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)
