Location: FCU_ExcitationContraction_Coupling @ 8bb03dd23f61 / python / FCU_EC_coupling.py

Author:
Shelley Fong <s.fong@auckland.ac.nz>
Date:
2022-02-18 14:17:34+13:00
Desc:
Updating way cellml is written
Permanent Source URI:
http://models.cellml.org/workspace/7a8/rawfile/8bb03dd23f61cba2ac9583321a845c02ff6f4155/python/FCU_EC_coupling.py

# ----------------------------------------------------------------
# This python script is the export of FCU_EC_coupling.cellml made on 4 Nov 2021
# ----------------------------------------------------------------
# Size of variable arrays:
sizeAlgebraic = 199
sizeStates = 49
sizeConstants = 143
from math import *
from numpy import *
import matplotlib.pyplot as plt
import time

def createLegends():
    legend_states = [""] * sizeStates
    legend_rates = [""] * sizeStates
    legend_algebraic = [""] * sizeAlgebraic
    legend_voi = ""
    legend_constants = [""] * sizeConstants
    legend_constants[0] = "kappa_PLBph1 in component BG_parameters (fmol_per_sec)"
    legend_constants[1] = "kappa_PLBph2 in component BG_parameters (fmol_per_sec)"
    legend_constants[2] = "kappa_PLBd1 in component BG_parameters (fmol_per_sec)"
    legend_constants[3] = "kappa_PLBd2 in component BG_parameters (fmol_per_sec)"
    legend_constants[4] = "kappa_Inh in component BG_parameters (fmol_per_sec)"
    legend_constants[5] = "kappa_SERCA_R1_2 in component BG_parameters (fmol_per_sec)"
    legend_constants[6] = "kappa_SERCA_R2_4 in component BG_parameters (fmol_per_sec)"
    legend_constants[7] = "kappa_SERCA_R2_2a in component BG_parameters (fmol_per_sec)"
    legend_constants[8] = "kappa_SERCA_R4_5 in component BG_parameters (fmol_per_sec)"
    legend_constants[9] = "kappa_SERCA_R5_6 in component BG_parameters (fmol_per_sec)"
    legend_constants[10] = "kappa_SERCA_R6_8 in component BG_parameters (fmol_per_sec)"
    legend_constants[11] = "kappa_SERCA_R8_9 in component BG_parameters (fmol_per_sec)"
    legend_constants[12] = "kappa_SERCA_R9_10 in component BG_parameters (fmol_per_sec)"
    legend_constants[13] = "kappa_SERCA_R10_1 in component BG_parameters (fmol_per_sec)"
    legend_constants[14] = "kappa_RyR in component BG_parameters (fmol_per_sec)"
    legend_constants[15] = "kappa_OC in component BG_parameters (fmol_per_sec)"
    legend_constants[16] = "kappa_CCI in component BG_parameters (fmol_per_sec)"
    legend_constants[17] = "kappa_CII in component BG_parameters (fmol_per_sec)"
    legend_constants[18] = "kappa_IO in component BG_parameters (fmol_per_sec)"
    legend_constants[19] = "kappa_Ca1 in component BG_parameters (fmol_per_sec)"
    legend_constants[20] = "kappa_Ca2 in component BG_parameters (fmol_per_sec)"
    legend_constants[21] = "kappa_K1 in component BG_parameters (fmol_per_sec)"
    legend_constants[22] = "kappa_K2 in component BG_parameters (fmol_per_sec)"
    legend_constants[23] = "kappa_d000 in component BG_parameters (fmol_per_sec)"
    legend_constants[24] = "kappa_d010 in component BG_parameters (fmol_per_sec)"
    legend_constants[25] = "kappa_d020 in component BG_parameters (fmol_per_sec)"
    legend_constants[26] = "kappa_d001 in component BG_parameters (fmol_per_sec)"
    legend_constants[27] = "kappa_d011 in component BG_parameters (fmol_per_sec)"
    legend_constants[28] = "kappa_d021 in component BG_parameters (fmol_per_sec)"
    legend_constants[29] = "kappa_f1_000 in component BG_parameters (fmol_per_sec)"
    legend_constants[30] = "kappa_f1_100 in component BG_parameters (fmol_per_sec)"
    legend_constants[31] = "kappa_f1_001 in component BG_parameters (fmol_per_sec)"
    legend_constants[32] = "kappa_f1_101 in component BG_parameters (fmol_per_sec)"
    legend_constants[33] = "kappa_f2_000 in component BG_parameters (fmol_per_sec)"
    legend_constants[34] = "kappa_f2_100 in component BG_parameters (fmol_per_sec)"
    legend_constants[35] = "kappa_f2_001 in component BG_parameters (fmol_per_sec)"
    legend_constants[36] = "kappa_f2_101 in component BG_parameters (fmol_per_sec)"
    legend_constants[37] = "kappa_f3_010 in component BG_parameters (fmol_per_sec)"
    legend_constants[38] = "kappa_f3_110 in component BG_parameters (fmol_per_sec)"
    legend_constants[39] = "kappa_f3_011 in component BG_parameters (fmol_per_sec)"
    legend_constants[40] = "kappa_f3_111 in component BG_parameters (fmol_per_sec)"
    legend_constants[41] = "kappa_fCa000 in component BG_parameters (fmol_per_sec)"
    legend_constants[42] = "kappa_fCa100 in component BG_parameters (fmol_per_sec)"
    legend_constants[43] = "kappa_fCa010 in component BG_parameters (fmol_per_sec)"
    legend_constants[44] = "kappa_fCa110 in component BG_parameters (fmol_per_sec)"
    legend_constants[45] = "kappa_fCa020 in component BG_parameters (fmol_per_sec)"
    legend_constants[46] = "kappa_fCa120 in component BG_parameters (fmol_per_sec)"
    legend_constants[47] = "kappa_R_TRPNCa in component BG_parameters (fmol_per_sec)"
    legend_constants[48] = "K_PLB in component BG_parameters (per_fmol)"
    legend_constants[49] = "K_PKACI in component BG_parameters (per_fmol)"
    legend_constants[50] = "K_PLB_PKACI in component BG_parameters (per_fmol)"
    legend_constants[51] = "K_PP1 in component BG_parameters (per_fmol)"
    legend_constants[52] = "K_PLBp_PP1 in component BG_parameters (per_fmol)"
    legend_constants[53] = "K_PLBp in component BG_parameters (per_fmol)"
    legend_constants[54] = "K_Ip in component BG_parameters (per_fmol)"
    legend_constants[55] = "K_Ip_PP1 in component BG_parameters (per_fmol)"
    legend_constants[56] = "K_P1_SERCA in component BG_parameters (per_fmol)"
    legend_constants[57] = "K_P2_SERCA in component BG_parameters (per_fmol)"
    legend_constants[58] = "K_P2a_SERCA in component BG_parameters (per_fmol)"
    legend_constants[59] = "K_P4_SERCA in component BG_parameters (per_fmol)"
    legend_constants[60] = "K_P5_SERCA in component BG_parameters (per_fmol)"
    legend_constants[61] = "K_P6_SERCA in component BG_parameters (per_fmol)"
    legend_constants[62] = "K_P8_SERCA in component BG_parameters (per_fmol)"
    legend_constants[63] = "K_P9_SERCA in component BG_parameters (per_fmol)"
    legend_constants[64] = "K_P10_SERCA in component BG_parameters (per_fmol)"
    legend_constants[65] = "K_H in component BG_parameters (per_fmol)"
    legend_constants[66] = "K_Cai in component BG_parameters (per_fmol)"
    legend_constants[67] = "K_Ca_SR in component BG_parameters (per_fmol)"
    legend_constants[68] = "K_MgATP in component BG_parameters (per_fmol)"
    legend_constants[69] = "K_MgADP in component BG_parameters (per_fmol)"
    legend_constants[70] = "K_P in component BG_parameters (per_fmol)"
    legend_constants[71] = "K_C_RyR in component BG_parameters (per_fmol)"
    legend_constants[72] = "K_CI_RyR in component BG_parameters (per_fmol)"
    legend_constants[73] = "K_I_RyR in component BG_parameters (per_fmol)"
    legend_constants[74] = "K_O_RyR in component BG_parameters (per_fmol)"
    legend_constants[75] = "K_Cao in component BG_parameters (per_fmol)"
    legend_constants[76] = "K_Ki in component BG_parameters (per_fmol)"
    legend_constants[77] = "K_Ko in component BG_parameters (per_fmol)"
    legend_constants[78] = "K_000_LCC in component BG_parameters (per_fmol)"
    legend_constants[79] = "K_010_LCC in component BG_parameters (per_fmol)"
    legend_constants[80] = "K_020_LCC in component BG_parameters (per_fmol)"
    legend_constants[81] = "K_100_LCC in component BG_parameters (per_fmol)"
    legend_constants[82] = "K_110_LCC in component BG_parameters (per_fmol)"
    legend_constants[83] = "K_120_LCC in component BG_parameters (per_fmol)"
    legend_constants[84] = "K_001_LCC in component BG_parameters (per_fmol)"
    legend_constants[85] = "K_011_LCC in component BG_parameters (per_fmol)"
    legend_constants[86] = "K_021_LCC in component BG_parameters (per_fmol)"
    legend_constants[87] = "K_101_LCC in component BG_parameters (per_fmol)"
    legend_constants[88] = "K_111_LCC in component BG_parameters (per_fmol)"
    legend_constants[89] = "K_121_LCC in component BG_parameters (per_fmol)"
    legend_constants[90] = "K_TRPN in component BG_parameters (per_fmol)"
    legend_constants[91] = "K_Ca_TRPN in component BG_parameters (per_fmol)"
    legend_voi = "time in component environment (second)"
    legend_constants[92] = "vol_myo in component environment (pL)"
    legend_constants[93] = "freq in component environment (dimensionless)"
    legend_constants[94] = "C_m in component environment (fF)"
    legend_states[0] = "q_membrane in component environment (fC)"
    legend_algebraic[35] = "V_m in component environment (volt)"
    legend_algebraic[63] = "v_RyR in component RyR (fmol_per_sec)"
    legend_constants[95] = "F in component constants (C_per_mol)"
    legend_constants[96] = "q_PLB_init in component environment (fmol)"
    legend_constants[97] = "q_PKACI_init in component environment (fmol)"
    legend_constants[98] = "q_PLB_PKACI_init in component environment (fmol)"
    legend_constants[99] = "q_PP1_init in component environment (fmol)"
    legend_constants[100] = "q_PLBp_PP1_init in compone.nt environment (fmol)"
    legend_constants[101] = "q_PLBp_init in component environment (fmol)"
    legend_constants[102] = "q_Ip_init in component environment (fmol)"
    legend_constants[103] = "q_Ip_PP1_init in component environment (fmol)"
    legend_constants[104] = "q_Cai_init in component environment (fmol)"
    legend_constants[105] = "q_Ca_SR_init in component environment (fmol)"
    legend_constants[106] = "q_H_init in component environment (fmol)"
    legend_constants[107] = "q_P_init in component environment (fmol)"
    legend_constants[108] = "q_MgADP_init in component environment (fmol)"
    legend_constants[109] = "q_MgATP_init in component environment (fmol)"
    legend_constants[110] = "q_Cao_init in component environment (fmol)"
    legend_constants[111] = "q_Ki_init in component environment (fmol)"
    legend_constants[112] = "q_Ko_init in component environment (fmol)"
    legend_constants[113] = "q_TRPN_init in component environment (fmol)"
    legend_constants[114] = "q_Ca_TRPN_init in component environment (fmol)"
    legend_algebraic[0] = "q_PLB in component environment (fmol)"
    legend_algebraic[8] = "q_PKACI in component environment (fmol)"
    legend_algebraic[10] = "q_PLB_PKACI in component environment (fmol)"
    legend_algebraic[12] = "q_PP1 in component environment (fmol)"
    legend_algebraic[16] = "q_PLBp_PP1 in component environment (fmol)"
    legend_algebraic[20] = "q_PLBp in component environment (fmol)"
    legend_algebraic[25] = "q_Ip in component environment (fmol)"
    legend_algebraic[32] = "q_Ip_PP1 in component environment (fmol)"
    legend_algebraic[1] = "q_H in component environment (fmol)"
    legend_algebraic[9] = "q_Cai in component environment (fmol)"
    legend_algebraic[13] = "q_MgATP in component environment (fmol)"
    legend_algebraic[17] = "q_MgADP in component environment (fmol)"
    legend_algebraic[21] = "q_P in component environment (fmol)"
    legend_algebraic[11] = "q_Ca_SR in component environment (fmol)"
    legend_algebraic[15] = "q_Cao in component environment (fmol)"
    legend_algebraic[19] = "q_Ki in component environment (fmol)"
    legend_algebraic[24] = "q_Ko in component environment (fmol)"
    legend_algebraic[18] = "q_TRPN in component environment (fmol)"
    legend_algebraic[22] = "q_Ca_TRPN in component environment (fmol)"
    legend_states[1] = "q_PLB in component PLB (fmol)"
    legend_states[2] = "q_PKACI in component PLB (fmol)"
    legend_states[3] = "q_PLB_PKACI in component PLB (fmol)"
    legend_states[4] = "q_PP1 in component PLB (fmol)"
    legend_states[5] = "q_PLBp_PP1 in component PLB (fmol)"
    legend_states[6] = "q_PLBp in component PLB (fmol)"
    legend_states[7] = "q_Ip in component PLB (fmol)"
    legend_states[8] = "q_Ip_PP1 in component PLB (fmol)"
    legend_states[9] = "q_H in component SERCA (fmol)"
    legend_states[10] = "q_Cai in component SERCA (fmol)"
    legend_states[11] = "q_Ca_SR in component SERCA (fmol)"
    legend_states[12] = "q_MgATP in component SERCA (fmol)"
    legend_states[13] = "q_MgADP in component SERCA (fmol)"
    legend_states[14] = "q_P in component SERCA (fmol)"
    legend_algebraic[26] = "q_SERCA_Ca_complexes in component SERCA (fmol)"
    legend_states[15] = "q_Ca_SR in component RyR (fmol)"
    legend_states[16] = "q_Cai in component RyR (fmol)"
    legend_algebraic[33] = "q_Ca_gate_complexes in component RyR_gate (fmol)"
    legend_states[17] = "q_Cai in component LCC (fmol)"
    legend_states[18] = "q_Cao in component LCC (fmol)"
    legend_states[19] = "q_Ki in component LCC (fmol)"
    legend_states[20] = "q_Ko in component LCC (fmol)"
    legend_algebraic[164] = "v_Ca_fgate in component LCC (fmol_per_sec)"
    legend_algebraic[38] = "q_Ca_fgate in component LCC_gate (fmol)"
    legend_states[21] = "q_Cai in component TRPN (fmol)"
    legend_states[22] = "q_TRPN in component TRPN (fmol)"
    legend_states[23] = "q_Ca_TRPN in component TRPN (fmol)"
    legend_algebraic[30] = "I_pulse in component environment (fA)"
    legend_constants[115] = "pulse_start in component environment (second)"
    legend_constants[116] = "pulse_end in component environment (second)"
    legend_constants[117] = "pulseMag in component environment (fA)"
    legend_constants[118] = "pulseHolding in component environment (fA)"
    legend_algebraic[67] = "v_E_RyR in component environment (fA)"
    legend_algebraic[197] = "v_E_LCC in component LCC (fA)"
    legend_algebraic[198] = "sum_I in component environment (fA)"
    legend_constants[119] = "zCa in component ion_valences (dimensionless)"
    legend_constants[120] = "zK in component ion_valences (dimensionless)"
    legend_algebraic[44] = "Ca_T in component environment (fmol)"
    legend_algebraic[27] = "PLB_T in component environment (fmol)"
    legend_algebraic[14] = "PKACI_T in component environment (fmol)"
    legend_algebraic[39] = "Ip_T in component environment (fmol)"
    legend_algebraic[23] = "K_i_T in component environment (fmol)"
    legend_algebraic[28] = "K_o_T in component environment (fmol)"
    legend_constants[121] = "R in component constants (J_per_K_per_mol)"
    legend_constants[122] = "T in component constants (kelvin)"
    legend_constants[123] = "zNa in component ion_valences (dimensionless)"
    legend_constants[124] = "zCl in component ion_valences (dimensionless)"
    legend_algebraic[68] = "vPLBph1 in component PLB (fmol_per_sec)"
    legend_algebraic[72] = "vPLBph2 in component PLB (fmol_per_sec)"
    legend_algebraic[76] = "vPLBd1 in component PLB (fmol_per_sec)"
    legend_algebraic[80] = "vPLBd2 in component PLB (fmol_per_sec)"
    legend_algebraic[84] = "vInh in component PLB (fmol_per_sec)"
    legend_algebraic[40] = "mu_PLB in component PLB (J_per_mol)"
    legend_algebraic[45] = "mu_PKACI in component PLB (J_per_mol)"
    legend_algebraic[49] = "mu_PLB_PKACI in component PLB (J_per_mol)"
    legend_algebraic[52] = "mu_PP1 in component PLB (J_per_mol)"
    legend_algebraic[55] = "mu_PLBp_PP1 in component PLB (J_per_mol)"
    legend_algebraic[58] = "mu_PLBp in component PLB (J_per_mol)"
    legend_algebraic[61] = "mu_Ip in component PLB (J_per_mol)"
    legend_algebraic[64] = "mu_Ip_PP1 in component PLB (J_per_mol)"
    legend_constants[125] = "n_Cai in component SERCA (dimensionless)"
    legend_constants[126] = "n_Ca_SR in component SERCA (dimensionless)"
    legend_constants[127] = "n_H in component SERCA (dimensionless)"
    legend_states[24] = "q_P1_SERCA in component SERCA (fmol)"
    legend_states[25] = "q_P2_SERCA in component SERCA (fmol)"
    legend_states[26] = "q_P2a_SERCA in component SERCA (fmol)"
    legend_states[27] = "q_P4_SERCA in component SERCA (fmol)"
    legend_states[28] = "q_P5_SERCA in component SERCA (fmol)"
    legend_states[29] = "q_P10_SERCA in component SERCA (fmol)"
    legend_states[30] = "q_P6_SERCA in component SERCA (fmol)"
    legend_states[31] = "q_P8_SERCA in component SERCA (fmol)"
    legend_states[32] = "q_P9_SERCA in component SERCA (fmol)"
    legend_algebraic[2] = "c_Cai in component SERCA (mM)"
    legend_algebraic[3] = "c_Ca_SR in component SERCA (mM)"
    legend_algebraic[4] = "c_H in component SERCA (mM)"
    legend_algebraic[5] = "c_MgADP in component SERCA (mM)"
    legend_algebraic[6] = "c_MgATP in component SERCA (mM)"
    legend_algebraic[7] = "c_P in component SERCA (mM)"
    legend_constants[128] = "vol_i in component SERCA (pL)"
    legend_constants[141] = "vol_sr in component SERCA (pL)"
    legend_constants[142] = "vol_isr in component SERCA (pL)"
    legend_algebraic[29] = "mu_Cai in component SERCA (J_per_mol)"
    legend_algebraic[139] = "v_Cai in component SERCA (fmol_per_sec)"
    legend_algebraic[34] = "mu_Ca_SR in component SERCA (J_per_mol)"
    legend_algebraic[146] = "v_Ca_SR in component SERCA (fmol_per_sec)"
    legend_algebraic[41] = "mu_H in component SERCA (J_per_mol)"
    legend_algebraic[156] = "v_H in component SERCA (fmol_per_sec)"
    legend_algebraic[46] = "mu_MgADP in component SERCA (J_per_mol)"
    legend_algebraic[133] = "v_MgADP in component SERCA (fmol_per_sec)"
    legend_algebraic[50] = "mu_MgATP in component SERCA (J_per_mol)"
    legend_algebraic[129] = "v_MgATP in component SERCA (fmol_per_sec)"
    legend_algebraic[56] = "mu_P1 in component SERCA (J_per_mol)"
    legend_algebraic[130] = "v_P1 in component SERCA (fmol_per_sec)"
    legend_algebraic[53] = "mu_P in component SERCA (J_per_mol)"
    legend_algebraic[125] = "v_P in component SERCA (fmol_per_sec)"
    legend_algebraic[59] = "mu_P2 in component SERCA (J_per_mol)"
    legend_algebraic[140] = "v_P2 in component SERCA (fmol_per_sec)"
    legend_algebraic[62] = "mu_P2a in component SERCA (J_per_mol)"
    legend_algebraic[136] = "v_P2a in component SERCA (fmol_per_sec)"
    legend_algebraic[65] = "mu_P4 in component SERCA (J_per_mol)"
    legend_algebraic[143] = "v_P4 in component SERCA (fmol_per_sec)"
    legend_algebraic[69] = "mu_P5 in component SERCA (J_per_mol)"
    legend_algebraic[144] = "v_P5 in component SERCA (fmol_per_sec)"
    legend_algebraic[77] = "mu_P6 in component SERCA (J_per_mol)"
    legend_algebraic[147] = "v_P6 in component SERCA (fmol_per_sec)"
    legend_algebraic[81] = "mu_P8 in component SERCA (J_per_mol)"
    legend_algebraic[151] = "v_P8 in component SERCA (fmol_per_sec)"
    legend_algebraic[85] = "mu_P9 in component SERCA (J_per_mol)"
    legend_algebraic[157] = "v_P9 in component SERCA (fmol_per_sec)"
    legend_algebraic[73] = "mu_P10 in component SERCA (J_per_mol)"
    legend_algebraic[158] = "v_P10 in component SERCA (fmol_per_sec)"
    legend_algebraic[91] = "Af_R1_2 in component SERCA (J_per_mol)"
    legend_algebraic[93] = "Ar_R1_2 in component SERCA (J_per_mol)"
    legend_algebraic[126] = "v_SERCA_R1_2 in component SERCA (fmol_per_sec)"
    legend_algebraic[95] = "Af_R5_6 in component SERCA (J_per_mol)"
    legend_algebraic[97] = "Ar_R5_6 in component SERCA (J_per_mol)"
    legend_algebraic[131] = "v_SERCA_R5_6 in component SERCA (fmol_per_sec)"
    legend_algebraic[99] = "Af_R2_2a in component SERCA (J_per_mol)"
    legend_algebraic[101] = "Ar_R2_2a in component SERCA (J_per_mol)"
    legend_algebraic[134] = "v_SERCA_R2_2a in component SERCA (fmol_per_sec)"
    legend_algebraic[103] = "Af_R2_4 in component SERCA (J_per_mol)"
    legend_algebraic[105] = "Ar_R2_4 in component SERCA (J_per_mol)"
    legend_algebraic[137] = "v_SERCA_R2_4 in component SERCA (fmol_per_sec)"
    legend_algebraic[107] = "Af_R4_5 in component SERCA (J_per_mol)"
    legend_algebraic[109] = "Ar_R4_5 in component SERCA (J_per_mol)"
    legend_algebraic[141] = "v_SERCA_R4_5 in component SERCA (fmol_per_sec)"
    legend_algebraic[111] = "Af_R6_8 in component SERCA (J_per_mol)"
    legend_algebraic[113] = "Ar_R6_8 in component SERCA (J_per_mol)"
    legend_algebraic[145] = "v_SERCA_R6_8 in component SERCA (fmol_per_sec)"
    legend_algebraic[115] = "Af_R8_9 in component SERCA (J_per_mol)"
    legend_algebraic[117] = "Ar_R8_9 in component SERCA (J_per_mol)"
    legend_algebraic[148] = "v_SERCA_R8_9 in component SERCA (fmol_per_sec)"
    legend_algebraic[119] = "Af_R9_10 in component SERCA (J_per_mol)"
    legend_algebraic[121] = "Ar_R9_10 in component SERCA (J_per_mol)"
    legend_algebraic[152] = "v_SERCA_R9_10 in component SERCA (fmol_per_sec)"
    legend_algebraic[87] = "Af_R10_1 in component SERCA (J_per_mol)"
    legend_algebraic[89] = "Ar_R10_1 in component SERCA (J_per_mol)"
    legend_algebraic[123] = "v_SERCA_R10_1 in component SERCA (fmol_per_sec)"
    legend_constants[129] = "nCa_1 in component RyR (dimensionless)"
    legend_constants[130] = "nCa_2 in component RyR (dimensionless)"
    legend_algebraic[42] = "mu_Ca_SR in component RyR (J_per_mol)"
    legend_algebraic[47] = "mu_Cai in component RyR (J_per_mol)"
    legend_algebraic[82] = "v_gate_Cai in component RyR (fmol_per_sec)"
    legend_algebraic[60] = "mu_O_RyR in component RyR_gate (J_per_mol)"
    legend_algebraic[66] = "v_OC in component RyR_gate (fmol_per_sec)"
    legend_algebraic[70] = "v_CCI in component RyR_gate (fmol_per_sec)"
    legend_algebraic[74] = "v_CII in component RyR_gate (fmol_per_sec)"
    legend_algebraic[78] = "v_IO in component RyR_gate (fmol_per_sec)"
    legend_states[33] = "q_C_RyR in component RyR_gate (fmol)"
    legend_states[34] = "q_CI_RyR in component RyR_gate (fmol)"
    legend_states[35] = "q_I_RyR in component RyR_gate (fmol)"
    legend_states[36] = "q_O_RyR in component RyR_gate (fmol)"
    legend_algebraic[51] = "mu_C_RyR in component RyR_gate (J_per_mol)"
    legend_algebraic[54] = "mu_CI_RyR in component RyR_gate (J_per_mol)"
    legend_algebraic[57] = "mu_I_RyR in component RyR_gate (J_per_mol)"
    legend_constants[140] = "RT in component LCC (J_per_mol)"
    legend_algebraic[71] = "mu_E in component LCC (J_per_C)"
    legend_algebraic[88] = "mu_mem_Ca in component LCC (J_per_mol)"
    legend_algebraic[90] = "mu_mem_K in component LCC (J_per_mol)"
    legend_algebraic[166] = "v_CaGHK_i in component LCC (fmol_per_sec)"
    legend_algebraic[149] = "v_CaGHK_o in component LCC (fmol_per_sec)"
    legend_algebraic[153] = "v_KGHK_i in component LCC (fmol_per_sec)"
    legend_algebraic[154] = "v_KGHK_o in component LCC (fmol_per_sec)"
    legend_algebraic[75] = "mu_Cai in component LCC (J_per_mol)"
    legend_algebraic[79] = "mu_Cao in component LCC (J_per_mol)"
    legend_algebraic[83] = "mu_Ki in component LCC (J_per_mol)"
    legend_algebraic[86] = "mu_Ko in component LCC (J_per_mol)"
    legend_algebraic[110] = "mu_S111 in component LCC_gate (J_per_mol)"
    legend_algebraic[127] = "mu_S121 in component LCC_gate (J_per_mol)"
    legend_algebraic[112] = "Af_Ca1 in component LCC (J_per_mol)"
    legend_algebraic[114] = "Ar_Ca1 in component LCC (J_per_mol)"
    legend_algebraic[120] = "v_Ca1 in component LCC (fmol_per_sec)"
    legend_algebraic[142] = "v_Ca2 in component LCC (fmol_per_sec)"
    legend_algebraic[122] = "v_K1 in component LCC (fmol_per_sec)"
    legend_algebraic[150] = "v_K2 in component LCC (fmol_per_sec)"
    legend_algebraic[128] = "Af_Ca2 in component LCC (J_per_mol)"
    legend_algebraic[132] = "Ar_Ca2 in component LCC (J_per_mol)"
    legend_algebraic[116] = "Af_K1 in component LCC (J_per_mol)"
    legend_algebraic[118] = "Ar_K1 in component LCC (J_per_mol)"
    legend_algebraic[135] = "Af_K2 in component LCC (J_per_mol)"
    legend_algebraic[138] = "Ar_K2 in component LCC (J_per_mol)"
    legend_algebraic[155] = "v_fCa00 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[159] = "v_fCa01 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[160] = "v_fCa02 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[161] = "v_fCa10 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[162] = "v_fCa11 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[163] = "v_fCa12 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[194] = "v_gate in component LCC_gate (fmol_per_sec)"
    legend_algebraic[36] = "v_pulse in component LCC (fmol_per_sec)"
    legend_states[37] = "q_S000 in component LCC_gate (fmol)"
    legend_states[38] = "q_S010 in component LCC_gate (fmol)"
    legend_states[39] = "q_S020 in component LCC_gate (fmol)"
    legend_states[40] = "q_S100 in component LCC_gate (fmol)"
    legend_states[41] = "q_S110 in component LCC_gate (fmol)"
    legend_states[42] = "q_S120 in component LCC_gate (fmol)"
    legend_states[43] = "q_S001 in component LCC_gate (fmol)"
    legend_states[44] = "q_S011 in component LCC_gate (fmol)"
    legend_states[45] = "q_S021 in component LCC_gate (fmol)"
    legend_states[46] = "q_S101 in component LCC_gate (fmol)"
    legend_states[47] = "q_S111 in component LCC_gate (fmol)"
    legend_states[48] = "q_S121 in component LCC_gate (fmol)"
    legend_algebraic[92] = "mu_S000 in component LCC_gate (J_per_mol)"
    legend_algebraic[96] = "mu_S010 in component LCC_gate (J_per_mol)"
    legend_algebraic[100] = "mu_S020 in component LCC_gate (J_per_mol)"
    legend_algebraic[104] = "mu_S100 in component LCC_gate (J_per_mol)"
    legend_algebraic[108] = "mu_S110 in component LCC_gate (J_per_mol)"
    legend_algebraic[124] = "mu_S120 in component LCC_gate (J_per_mol)"
    legend_algebraic[94] = "mu_S001 in component LCC_gate (J_per_mol)"
    legend_algebraic[98] = "mu_S011 in component LCC_gate (J_per_mol)"
    legend_algebraic[102] = "mu_S021 in component LCC_gate (J_per_mol)"
    legend_algebraic[106] = "mu_S101 in component LCC_gate (J_per_mol)"
    legend_algebraic[172] = "v_f1_000 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[173] = "v_f2_000 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[174] = "v_f3_010 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[176] = "v_f1_100 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[178] = "v_f2_100 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[180] = "v_f3_110 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[165] = "v_d000 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[167] = "v_d010 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[168] = "v_d020 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[169] = "v_d001 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[170] = "v_d011 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[171] = "v_d021 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[182] = "v_f1_001 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[185] = "v_f2_001 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[186] = "v_f3_011 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[188] = "v_f1_101 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[190] = "v_f2_101 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[192] = "v_f3_111 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[175] = "v_S000 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[179] = "v_S010 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[177] = "v_S020 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[181] = "v_S100 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[183] = "v_S110 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[184] = "v_S120 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[187] = "v_S001 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[191] = "v_S011 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[189] = "v_S021 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[193] = "v_S101 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[195] = "v_S111 in component LCC_gate (fmol_per_sec)"
    legend_algebraic[196] = "v_S121 in component LCC_gate (fmol_per_sec)"
    legend_constants[131] = "z_rCa in component LCC_gate (dimensionless)"
    legend_constants[132] = "z_fd in component LCC_gate (dimensionless)"
    legend_constants[133] = "z_ff1 in component LCC_gate (dimensionless)"
    legend_constants[134] = "z_ff2 in component LCC_gate (dimensionless)"
    legend_constants[135] = "z_ff3 in component LCC_gate (dimensionless)"
    legend_constants[136] = "z_rd in component LCC_gate (dimensionless)"
    legend_constants[137] = "z_rf1 in component LCC_gate (dimensionless)"
    legend_constants[138] = "z_rf2 in component LCC_gate (dimensionless)"
    legend_constants[139] = "z_rf3 in component LCC_gate (dimensionless)"
    legend_algebraic[31] = "mu_Cai in component TRPN (J_per_mol)"
    legend_algebraic[37] = "mu_TRPN in component TRPN (J_per_mol)"
    legend_algebraic[43] = "mu_Ca_TRPN in component TRPN (J_per_mol)"
    legend_algebraic[48] = "v_R_TRPNCa in component TRPN (fmol_per_sec)"
    legend_rates[0] = "d/dt q_membrane in component environment (fC)"
    legend_rates[1] = "d/dt q_PLB in component PLB (fmol)"
    legend_rates[2] = "d/dt q_PKACI in component PLB (fmol)"
    legend_rates[3] = "d/dt q_PLB_PKACI in component PLB (fmol)"
    legend_rates[4] = "d/dt q_PP1 in component PLB (fmol)"
    legend_rates[5] = "d/dt q_PLBp_PP1 in component PLB (fmol)"
    legend_rates[6] = "d/dt q_PLBp in component PLB (fmol)"
    legend_rates[7] = "d/dt q_Ip in component PLB (fmol)"
    legend_rates[8] = "d/dt q_Ip_PP1 in component PLB (fmol)"
    legend_rates[10] = "d/dt q_Cai in component SERCA (fmol)"
    legend_rates[11] = "d/dt q_Ca_SR in component SERCA (fmol)"
    legend_rates[9] = "d/dt q_H in component SERCA (fmol)"
    legend_rates[13] = "d/dt q_MgADP in component SERCA (fmol)"
    legend_rates[12] = "d/dt q_MgATP in component SERCA (fmol)"
    legend_rates[14] = "d/dt q_P in component SERCA (fmol)"
    legend_rates[24] = "d/dt q_P1_SERCA in component SERCA (fmol)"
    legend_rates[25] = "d/dt q_P2_SERCA in component SERCA (fmol)"
    legend_rates[26] = "d/dt q_P2a_SERCA in component SERCA (fmol)"
    legend_rates[27] = "d/dt q_P4_SERCA in component SERCA (fmol)"
    legend_rates[28] = "d/dt q_P5_SERCA in component SERCA (fmol)"
    legend_rates[29] = "d/dt q_P10_SERCA in component SERCA (fmol)"
    legend_rates[30] = "d/dt q_P6_SERCA in component SERCA (fmol)"
    legend_rates[31] = "d/dt q_P8_SERCA in component SERCA (fmol)"
    legend_rates[32] = "d/dt q_P9_SERCA in component SERCA (fmol)"
    legend_rates[15] = "d/dt q_Ca_SR in component RyR (fmol)"
    legend_rates[16] = "d/dt q_Cai in component RyR (fmol)"
    legend_rates[36] = "d/dt q_O_RyR in component RyR_gate (fmol)"
    legend_rates[33] = "d/dt q_C_RyR in component RyR_gate (fmol)"
    legend_rates[34] = "d/dt q_CI_RyR in component RyR_gate (fmol)"
    legend_rates[35] = "d/dt q_I_RyR in component RyR_gate (fmol)"
    legend_rates[17] = "d/dt q_Cai in component LCC (fmol)"
    legend_rates[18] = "d/dt q_Cao in component LCC (fmol)"
    legend_rates[19] = "d/dt q_Ki in component LCC (fmol)"
    legend_rates[20] = "d/dt q_Ko in component LCC (fmol)"
    legend_rates[37] = "d/dt q_S000 in component LCC_gate (fmol)"
    legend_rates[43] = "d/dt q_S001 in component LCC_gate (fmol)"
    legend_rates[38] = "d/dt q_S010 in component LCC_gate (fmol)"
    legend_rates[44] = "d/dt q_S011 in component LCC_gate (fmol)"
    legend_rates[39] = "d/dt q_S020 in component LCC_gate (fmol)"
    legend_rates[45] = "d/dt q_S021 in component LCC_gate (fmol)"
    legend_rates[40] = "d/dt q_S100 in component LCC_gate (fmol)"
    legend_rates[46] = "d/dt q_S101 in component LCC_gate (fmol)"
    legend_rates[41] = "d/dt q_S110 in component LCC_gate (fmol)"
    legend_rates[47] = "d/dt q_S111 in component LCC_gate (fmol)"
    legend_rates[42] = "d/dt q_S120 in component LCC_gate (fmol)"
    legend_rates[48] = "d/dt q_S121 in component LCC_gate (fmol)"
    legend_rates[21] = "d/dt q_Cai in component TRPN (fmol)"
    legend_rates[22] = "d/dt q_TRPN in component TRPN (fmol)"
    legend_rates[23] = "d/dt q_Ca_TRPN in component TRPN (fmol)"
    return (legend_states, legend_algebraic, legend_voi, legend_constants)

def initConsts(Gmult,igs):
    hackMultiplier = 1e0
    constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
    constants[0] = 45.5263
    constants[1] = 6.55904
    constants[2] = 0.386674
    constants[3] = 1.21269
    constants[4] = 431.435
    constants[5] = 2.24963e-05
    constants[6] = 108329
    constants[7] = 2.77146e+06
    constants[8] = 108329
    constants[9] = 3442.07
    constants[10] = 3.07427e+06
    constants[11] = 1.0306e+07
    constants[12] = 1.0306e+07
    constants[13] = 0.00565005
    constants[14] = 3.58172e+06
    constants[15] = 8.83262
    constants[16] = 0.12618
    constants[17] = 883.262
    constants[18] = 73.6052
    constants[19] = 8.37314
    constants[20] = 15.6575
    constants[21] = 0.0929479
    constants[22] = 0.17381
    constants[23] = 942.411
    constants[24] = 23.2508
    constants[25] = 43.4783
    constants[26] = 0.175516
    constants[27] = 0.00433026
    constants[28] = 0.00809745
    constants[29] = 1.73616
    constants[30] = 8.62132
    constants[31] = 0.000323344
    constants[32] = 0.00160565
    constants[33] = 3.17968
    constants[34] = 15.7895
    constants[35] = 0.000592187
    constants[36] = 0.00294065
    constants[37] = 8932.14
    constants[38] = 44354.8
    constants[39] = 1.66353
    constants[40] = 8.2607
    constants[41] = 60096.5
    constants[42] = 298424
    constants[43] = 1482.67
    constants[44] = 7362.59
    constants[45] = 2772.56
    constants[46] = 13767.8
    constants[47] = 223.356
    constants[48] = 0.00235741
    constants[49] = 0.395191
    constants[50] = 0.638527
    constants[51] = 0.361989
    constants[52] = 0.203757
    constants[53] = 0.326014
    constants[54] = 5.41093
    constants[55] = 0.0673793
    constants[56] = 906.879
    constants[57] = 2605.18
    constants[58] = 271.664
    constants[59] = 6950.17
    constants[60] = 0.0214514
    constants[61] = 1228.84
    constants[62] = 76.0193
    constants[63] = 73.055
    constants[64] = 105.473
    constants[65] = 278.106
    constants[66] = 0.052177
    constants[67] = 0.052177
    constants[68] = 1072.81
    constants[69] = 1.34257e-05
    constants[70] = 0.0245735
    constants[71] = 115.191
    constants[72] = 1.15191
    constants[73] = 0.00197471
    constants[74] = 0.197471
    constants[75] = 0.052177
    constants[76] = 0.00369486
    constants[77] = 0.00369486
    constants[78] = 0.0150147
    constants[79] = 0.608585
    constants[80] = 0.325451
    constants[81] = 0.00302366
    constants[82] = 0.122556
    constants[83] = 0.0655392
    constants[84] = 80.6198
    constants[85] = 3267.72
    constants[86] = 1747.47
    constants[87] = 16.2352
    constants[88] = 658.051
    constants[89] = 351.905
    constants[90] = 7.25115
    constants[91] = 0.02603
    constants[92] = 34.4
    constants[93] = 500
    constants[94] = 1.381e5
    states[0] = -13039
    constants[95] = 96485
    constants[96] = 4.028E+00
    constants[97] = 2.234E-03
    constants[98] = 1e-18
    constants[99] = 3.382E-02
    constants[100] = 1e-18
    constants[101] = 1e-18
    constants[102] = 1.999E-03
    constants[103] = 1e-18
    constants[104] = 6.82e-3
    constants[105] = 0.641
    constants[106] = 0.004028
    constants[107] = 570
    constants[108] = 1.3794
    constants[109] = 3.8
    constants[110] = 6.84
    constants[111] = 5.51E+03
    constants[112] = 2.05E+02
    constants[113] = 2.57
    constants[114] = 1e-18
    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-16
    states[19] = 1e-16
    states[20] = 1e-16
    states[21] = 1e-16
    states[22] = 1e-16
    states[23] = 1e-16
    constants[115] = 4e-2
    constants[116] = 4.1e-2
    constants[117] = 1e6
    constants[118] = 0
    constants[119] = 2
    constants[120] = 1
    constants[121] = 8.31
    constants[122] = 310
    constants[123] = 1
    constants[124] = -1
    constants[125] = 2
    constants[126] = 2
    constants[127] = 2
    states[24] = 0.11111
    states[25] = 0.11111
    states[26] = 0.11111
    states[27] = 0.11111
    states[28] = 0.11111
    states[29] = 0.11111
    states[30] = 0.11111
    states[31] = 0.11111
    states[32] = 0.11111
    constants[128] = 34.0
    constants[129] = 1
    constants[130] = 2
    states[33] = 0.25e-3
    states[34] = 0.25e-3
    states[35] = 0.25e-3
    states[36] = 0.25e-3
    states[37] = 3.6989e-07
    states[38] = 7.3239e-05
    states[39] = 3.6989e-07
    states[40] = 3.7363e-09
    states[41] = 7.3979e-07
    states[42] = 3.7363e-09
    states[43] = 4.1099e-08
    states[44] = 8.1377e-06
    states[45] = 4.1099e-08
    states[46] = 4.1514e-10
    states[47] = 8.2199e-08
    states[48] = 4.1514e-10
    constants[131] = 2
    constants[132] = 2.1404
    constants[133] = -1.1495
    constants[134] = 0.72162
    constants[135] = 4.2933
    constants[136] = -2.1404
    constants[137] = 1.8993
    constants[138] = -0.52288
    constants[139] = 0
    constants[140] = constants[121]*constants[122]
    constants[141] = constants[128]*0.0350000
    constants[142] = constants[128]+constants[141]
    # hack increase constants specified by igs by a multiplier Gmult
    constants = [c*Gmult if i in igs else c for i, c in enumerate(constants)]
    states = [s*hackMultiplier for s in states]
    return (states, constants)

def computeRates(voi, states, constants):
    print(voi)
    rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
    algebraic[9] = constants[104]+states[10]+states[17]+states[21]+states[16]
    algebraic[31] = constants[121]*constants[122]*log(constants[66]*algebraic[9])
    algebraic[18] = constants[113]+states[22]
    algebraic[37] = constants[121]*constants[122]*log(constants[90]*algebraic[18])
    algebraic[22] = constants[114]+states[23]
    algebraic[43] = constants[121]*constants[122]*log(constants[91]*algebraic[22])
    algebraic[48] = constants[47]*(exp((algebraic[31]+algebraic[37])/(constants[121]*constants[122]))-exp(algebraic[43]/(constants[121]*constants[122])))
    rates[21] = -algebraic[48]
    rates[22] = -algebraic[48]
    rates[23] = algebraic[48]
    algebraic[35] = states[0]/constants[94]
    algebraic[11] = constants[105]+states[11]+states[15]
    algebraic[42] = constants[121]*constants[122]*log(constants[67]*algebraic[11])
    algebraic[47] = constants[121]*constants[122]*log(constants[66]*algebraic[9])
    algebraic[60] = constants[121]*constants[122]*log(constants[74]*states[36])
    algebraic[63] = constants[14]*exp(algebraic[60]/(constants[121]*constants[122]))*(exp((algebraic[42]+constants[119]*constants[95]*algebraic[35])/(constants[121]*constants[122]))-exp(algebraic[47]/(constants[121]*constants[122])))
    rates[15] = -algebraic[63]
    algebraic[51] = constants[121]*constants[122]*log(constants[71]*states[33])
    algebraic[66] = constants[15]*(exp(algebraic[60]/(constants[121]*constants[122]))-exp((algebraic[51]+constants[130]*algebraic[47])/(constants[121]*constants[122])))
    algebraic[54] = constants[121]*constants[122]*log(constants[72]*states[34])
    algebraic[70] = constants[16]*(exp((algebraic[51]+constants[129]*algebraic[47])/(constants[121]*constants[122]))-exp(algebraic[54]/(constants[121]*constants[122])))
    rates[33] = algebraic[66]-algebraic[70]
    algebraic[0] = constants[96]+states[1]
    algebraic[40] = constants[121]*constants[122]*log(constants[48]*algebraic[0])
    algebraic[8] = constants[97]+states[2]
    algebraic[45] = constants[121]*constants[122]*log(constants[49]*algebraic[8])
    algebraic[10] = constants[98]+states[3]
    algebraic[49] = constants[121]*constants[122]*log(constants[50]*algebraic[10])
    algebraic[68] = constants[0]*(exp((algebraic[40]+algebraic[45])/(constants[121]*constants[122]))-exp(algebraic[49]/(constants[121]*constants[122])))
    algebraic[20] = constants[101]+states[6]
    algebraic[58] = constants[121]*constants[122]*log(constants[53]*algebraic[20])
    algebraic[72] = constants[1]*(exp(algebraic[49]/(constants[121]*constants[122]))-exp((algebraic[58]+algebraic[45])/(constants[121]*constants[122])))
    rates[2] = algebraic[72]-algebraic[68]
    rates[3] = algebraic[68]-algebraic[72]
    algebraic[57] = constants[121]*constants[122]*log(constants[73]*states[35])
    algebraic[74] = constants[17]*(exp((algebraic[54]+constants[130]*algebraic[47])/(constants[121]*constants[122]))-exp(algebraic[57]/(constants[121]*constants[122])))
    rates[34] = algebraic[70]-algebraic[74]
    algebraic[12] = constants[99]+states[4]
    algebraic[52] = constants[121]*constants[122]*log(constants[51]*algebraic[12])
    algebraic[16] = constants[100]+states[5]
    algebraic[55] = constants[121]*constants[122]*log(constants[52]*algebraic[16])
    algebraic[76] = constants[2]*(exp((algebraic[58]+algebraic[52])/(constants[121]*constants[122]))-exp(algebraic[55]/(constants[121]*constants[122])))
    rates[6] = algebraic[72]-algebraic[76]
    algebraic[78] = constants[18]*(exp(algebraic[57]/(constants[121]*constants[122]))-exp((algebraic[60]+constants[129]*algebraic[47])/(constants[121]*constants[122])))
    rates[36] = algebraic[78]-algebraic[66]
    rates[35] = algebraic[74]-algebraic[78]
    algebraic[80] = constants[3]*(exp(algebraic[55]/(constants[121]*constants[122]))-exp((algebraic[40]+algebraic[52])/(constants[121]*constants[122])))
    rates[1] = algebraic[80]-algebraic[68]
    rates[5] = algebraic[76]-algebraic[80]
    algebraic[82] = ((constants[130]*algebraic[66]-constants[129]*algebraic[70])-constants[130]*algebraic[74])+constants[129]*algebraic[78]
    rates[16] = algebraic[63]+algebraic[82]
    algebraic[25] = constants[102]+states[7]
    algebraic[61] = constants[121]*constants[122]*log(constants[54]*algebraic[25])
    algebraic[32] = constants[103]+states[8]
    algebraic[64] = constants[121]*constants[122]*log(constants[55]*algebraic[32])
    algebraic[84] = constants[4]*(exp((algebraic[52]+algebraic[61])/(constants[121]*constants[122]))-exp(algebraic[64]/(constants[121]*constants[122])))
    rates[4] = (algebraic[80]-algebraic[76])-algebraic[84]
    rates[7] = -algebraic[84]
    rates[8] = algebraic[84]
    algebraic[73] = constants[121]*constants[122]*log(constants[64]*states[29])
    algebraic[87] = algebraic[73]
    algebraic[56] = constants[121]*constants[122]*log(constants[56]*states[24])
    algebraic[21] = constants[107]+states[14]
    algebraic[53] = constants[121]*constants[122]*log(constants[70]*algebraic[21])
    algebraic[89] = algebraic[56]+algebraic[53]
    algebraic[123] = constants[13]*(exp(algebraic[87]/(constants[121]*constants[122]))-exp(algebraic[89]/(constants[121]*constants[122])))
    algebraic[125] = algebraic[123]
    rates[14] = algebraic[125]
    algebraic[13] = constants[109]+states[12]
    algebraic[50] = constants[121]*constants[122]*log(constants[68]*algebraic[13])
    algebraic[91] = algebraic[56]+algebraic[50]
    algebraic[59] = constants[121]*constants[122]*log(constants[57]*states[25])
    algebraic[93] = algebraic[59]
    algebraic[126] = constants[5]*(exp(algebraic[91]/(constants[121]*constants[122]))-exp(algebraic[93]/(constants[121]*constants[122])))
    algebraic[129] = -algebraic[126]
    rates[12] = algebraic[129]
    algebraic[130] = algebraic[123]-algebraic[126]
    rates[24] = algebraic[130]
    algebraic[69] = constants[121]*constants[122]*log(constants[60]*states[28])
    algebraic[95] = algebraic[69]
    algebraic[17] = constants[108]+states[13]
    algebraic[46] = constants[121]*constants[122]*log(constants[69]*algebraic[17])
    algebraic[77] = constants[121]*constants[122]*log(constants[61]*states[30])
    algebraic[97] = algebraic[46]+algebraic[77]
    algebraic[131] = constants[9]*(exp(algebraic[95]/(constants[121]*constants[122]))-exp(algebraic[97]/(constants[121]*constants[122])))
    algebraic[133] = algebraic[131]
    rates[13] = algebraic[133]
    algebraic[1] = constants[106]+states[9]
    algebraic[41] = constants[121]*constants[122]*log(constants[65]*algebraic[1])
    algebraic[99] = algebraic[59]+algebraic[41]
    algebraic[62] = constants[121]*constants[122]*log(constants[58]*states[26])
    algebraic[101] = algebraic[62]
    algebraic[134] = constants[7]*(exp(algebraic[99]/(constants[121]*constants[122]))-exp(algebraic[101]/(constants[121]*constants[122])))
    algebraic[136] = algebraic[134]
    rates[26] = algebraic[136]
    algebraic[29] = constants[121]*constants[122]*log(constants[66]*algebraic[9])
    algebraic[103] = algebraic[59]+constants[125]*algebraic[29]
    algebraic[65] = constants[121]*constants[122]*log(constants[59]*states[27])
    algebraic[105] = algebraic[65]
    algebraic[137] = constants[6]*(exp(algebraic[103]/(constants[121]*constants[122]))-exp(algebraic[105]/(constants[121]*constants[122])))
    algebraic[139] = -constants[125]*algebraic[137]
    rates[10] = algebraic[139]
    algebraic[140] = algebraic[126]-algebraic[137]
    rates[25] = algebraic[140]
    algebraic[107] = algebraic[65]
    algebraic[109] = algebraic[69]+constants[127]*algebraic[41]
    algebraic[141] = constants[8]*(exp(algebraic[107]/(constants[121]*constants[122]))-exp(algebraic[109]/(constants[121]*constants[122])))
    algebraic[143] = algebraic[137]-algebraic[141]
    rates[27] = algebraic[143]
    algebraic[144] = algebraic[141]-algebraic[131]
    rates[28] = algebraic[144]
    algebraic[111] = algebraic[77]
    algebraic[34] = constants[121]*constants[122]*log(constants[67]*algebraic[11])
    algebraic[81] = constants[121]*constants[122]*log(constants[62]*states[31])
    algebraic[113] = algebraic[81]+constants[126]*algebraic[34]
    algebraic[145] = constants[10]*(exp(algebraic[111]/(constants[121]*constants[122]))-exp(algebraic[113]/(constants[121]*constants[122])))
    algebraic[146] = constants[126]*algebraic[145]
    rates[11] = algebraic[146]
    algebraic[147] = algebraic[131]-algebraic[145]
    rates[30] = algebraic[147]
    algebraic[71] = states[0]/constants[94]
    algebraic[88] = constants[119]*constants[95]*algebraic[71]
    algebraic[75] = constants[140]*log(constants[66]*algebraic[9])
    algebraic[110] = constants[140]*log(constants[88]*states[47])
    algebraic[112] = algebraic[75]+algebraic[88]+algebraic[110]
    algebraic[15] = constants[110]+states[18]
    algebraic[79] = constants[140]*log(constants[75]*algebraic[15])
    algebraic[114] = algebraic[79]+algebraic[110]
    algebraic[120] = custom_piecewise([equal(algebraic[88] , 0.000000), constants[19]*(exp(algebraic[112]/constants[140])-exp(algebraic[114]/constants[140])) , True, (((constants[19]*algebraic[88])/constants[140])/(exp(algebraic[88]/constants[140])-1.00000))*(exp(algebraic[112]/constants[140])-exp(algebraic[114]/constants[140]))])
    algebraic[127] = constants[140]*log(constants[89]*states[48])
    algebraic[128] = algebraic[75]+algebraic[88]+algebraic[127]
    algebraic[132] = algebraic[79]+algebraic[127]
    algebraic[142] = custom_piecewise([equal(algebraic[88] , 0.000000), constants[20]*(exp(algebraic[128]/constants[140])-exp(algebraic[132]/constants[140])) , True, (((constants[20]*algebraic[88])/constants[140])/(exp(algebraic[88]/constants[140])-1.00000))*(exp(algebraic[128]/constants[140])-exp(algebraic[132]/constants[140]))])
    algebraic[149] = algebraic[142]+algebraic[120]
    rates[18] = algebraic[149]
    algebraic[115] = algebraic[81]+constants[127]*algebraic[41]
    algebraic[85] = constants[121]*constants[122]*log(constants[63]*states[32])
    algebraic[117] = algebraic[85]
    algebraic[148] = constants[11]*(exp(algebraic[115]/(constants[121]*constants[122]))-exp(algebraic[117]/(constants[121]*constants[122])))
    algebraic[151] = algebraic[145]-algebraic[148]
    rates[31] = algebraic[151]
    algebraic[90] = constants[120]*constants[95]*algebraic[71]
    algebraic[19] = constants[111]+states[19]
    algebraic[83] = constants[140]*log(constants[76]*algebraic[19])
    algebraic[116] = algebraic[83]+algebraic[90]+algebraic[110]
    algebraic[24] = constants[112]+states[20]
    algebraic[86] = constants[140]*log(constants[77]*algebraic[24])
    algebraic[118] = algebraic[86]+algebraic[110]
    algebraic[122] = custom_piecewise([equal(algebraic[90] , 0.000000), constants[21]*(exp(algebraic[116]/constants[140])-exp(algebraic[118]/constants[140])) , True, (((constants[21]*algebraic[90])/constants[140])/(exp(algebraic[90]/constants[140])-1.00000))*(exp(algebraic[116]/constants[140])-exp(algebraic[118]/constants[140]))])
    algebraic[135] = algebraic[83]+algebraic[90]+algebraic[127]
    algebraic[138] = algebraic[86]+algebraic[127]
    algebraic[150] = custom_piecewise([equal(algebraic[90] , 0.000000), constants[22]*(exp(algebraic[135]/constants[140])-exp(algebraic[138]/constants[140])) , True, (((constants[22]*algebraic[90])/constants[140])/(exp(algebraic[90]/constants[140])-1.00000))*(exp(algebraic[135]/constants[140])-exp(algebraic[138]/constants[140]))])
    algebraic[153] = -algebraic[122]-algebraic[150]
    rates[19] = algebraic[153]
    algebraic[154] = algebraic[150]+algebraic[122]
    rates[20] = algebraic[154]
    algebraic[119] = algebraic[85]
    algebraic[121] = algebraic[41]+algebraic[73]
    algebraic[152] = constants[12]*(exp(algebraic[119]/(constants[121]*constants[122]))-exp(algebraic[121]/(constants[121]*constants[122])))
    algebraic[156] = ((constants[127]*algebraic[141]-algebraic[134])-constants[127]*algebraic[148])+algebraic[152]
    rates[9] = algebraic[156]
    algebraic[158] = algebraic[152]-algebraic[123]
    rates[29] = algebraic[158]
    algebraic[157] = algebraic[148]-algebraic[152]
    rates[32] = algebraic[157]
    algebraic[92] = constants[140]*log(constants[78]*states[37])
    algebraic[94] = constants[140]*log(constants[84]*states[43])
    algebraic[155] = constants[41]*(exp(algebraic[92]/constants[140])-exp((algebraic[94]+constants[131]*algebraic[75])/constants[140]))
    algebraic[96] = constants[140]*log(constants[79]*states[38])
    algebraic[98] = constants[140]*log(constants[85]*states[44])
    algebraic[159] = constants[43]*(exp(algebraic[96]/constants[140])-exp((algebraic[98]+constants[131]*algebraic[75])/constants[140]))
    algebraic[100] = constants[140]*log(constants[80]*states[39])
    algebraic[102] = constants[140]*log(constants[86]*states[45])
    algebraic[160] = constants[45]*(exp(algebraic[100]/constants[140])-exp((algebraic[102]+constants[131]*algebraic[75])/constants[140]))
    algebraic[104] = constants[140]*log(constants[81]*states[40])
    algebraic[106] = constants[140]*log(constants[87]*states[46])
    algebraic[161] = constants[42]*(exp(algebraic[104]/constants[140])-exp((algebraic[106]+constants[131]*algebraic[75])/constants[140]))
    algebraic[108] = constants[140]*log(constants[82]*states[41])
    algebraic[162] = constants[44]*(exp(algebraic[108]/constants[140])-exp((algebraic[110]+constants[131]*algebraic[75])/constants[140]))
    algebraic[124] = constants[140]*log(constants[83]*states[42])
    algebraic[163] = constants[46]*(exp(algebraic[124]/constants[140])-exp((algebraic[127]+constants[131]*algebraic[75])/constants[140]))
    algebraic[164] = constants[119]*(algebraic[161]+algebraic[162]+algebraic[160]+algebraic[159]+algebraic[155]+algebraic[163])
    algebraic[166] = (algebraic[164]-algebraic[120])-algebraic[142]
    rates[17] = algebraic[166]
    algebraic[172] = constants[29]*(exp((algebraic[92]+constants[133]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[96]+constants[137]*constants[95]*algebraic[71])/constants[140]))
    algebraic[173] = constants[33]*(exp((algebraic[92]+constants[134]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[100]+constants[138]*constants[95]*algebraic[71])/constants[140]))
    algebraic[165] = constants[23]*(exp((algebraic[92]+constants[132]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[104]+constants[136]*constants[95]*algebraic[71])/constants[140]))
    algebraic[175] = ((-algebraic[165]-algebraic[172])-algebraic[173])-algebraic[155]
    rates[37] = algebraic[175]
    algebraic[174] = constants[37]*(exp((algebraic[96]+constants[135]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[100]+constants[139]*constants[95]*algebraic[71])/constants[140]))
    algebraic[168] = constants[25]*(exp((algebraic[100]+constants[132]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[124]+constants[136]*constants[95]*algebraic[71])/constants[140]))
    algebraic[177] = (-algebraic[168]+algebraic[174]+algebraic[173])-algebraic[160]
    rates[39] = algebraic[177]
    algebraic[176] = constants[30]*(exp((algebraic[104]+constants[133]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[108]+constants[137]*constants[95]*algebraic[71])/constants[140]))
    algebraic[167] = constants[24]*(exp((algebraic[96]+constants[132]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[108]+constants[136]*constants[95]*algebraic[71])/constants[140]))
    algebraic[179] = ((-algebraic[167]-algebraic[174])+algebraic[176])-algebraic[159]
    rates[38] = algebraic[179]
    algebraic[178] = constants[34]*(exp((algebraic[104]+constants[134]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[124]+constants[138]*constants[95]*algebraic[71])/constants[140]))
    algebraic[181] = ((-algebraic[178]+algebraic[165])-algebraic[176])-algebraic[161]
    rates[40] = algebraic[181]
    algebraic[180] = constants[38]*(exp((algebraic[108]+constants[135]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[124]+constants[139]*constants[95]*algebraic[71])/constants[140]))
    algebraic[183] = ((algebraic[176]+algebraic[167])-algebraic[162])-algebraic[180]
    rates[41] = algebraic[183]
    algebraic[184] = ((algebraic[178]+algebraic[168])-algebraic[163])+algebraic[180]
    rates[42] = algebraic[184]
    algebraic[169] = constants[26]*(exp((algebraic[94]+constants[132]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[106]+constants[136]*constants[95]*algebraic[71])/constants[140]))
    algebraic[182] = constants[31]*(exp((algebraic[94]+constants[133]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[98]+constants[137]*constants[95]*algebraic[71])/constants[140]))
    algebraic[185] = constants[35]*(exp((algebraic[94]+constants[134]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[102]+constants[138]*constants[95]*algebraic[71])/constants[140]))
    algebraic[187] = ((-algebraic[169]-algebraic[182])-algebraic[185])+algebraic[155]
    rates[43] = algebraic[187]
    algebraic[171] = constants[28]*(exp((algebraic[102]+constants[132]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[127]+constants[136]*constants[95]*algebraic[71])/constants[140]))
    algebraic[186] = constants[39]*(exp((algebraic[98]+constants[135]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[102]+constants[139]*constants[95]*algebraic[71])/constants[140]))
    algebraic[189] = -algebraic[171]+algebraic[186]+algebraic[185]+algebraic[160]
    rates[45] = algebraic[189]
    algebraic[170] = constants[27]*(exp((algebraic[98]+constants[132]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[110]+constants[136]*constants[95]*algebraic[71])/constants[140]))
    algebraic[188] = constants[32]*(exp((algebraic[106]+constants[133]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[110]+constants[137]*constants[95]*algebraic[71])/constants[140]))
    algebraic[191] = (-algebraic[170]-algebraic[186])+algebraic[188]+algebraic[159]
    rates[44] = algebraic[191]
    algebraic[190] = constants[36]*(exp((algebraic[106]+constants[134]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[127]+constants[138]*constants[95]*algebraic[71])/constants[140]))
    algebraic[193] = ((-algebraic[190]+algebraic[169])-algebraic[188])+algebraic[161]
    rates[46] = algebraic[193]
    algebraic[192] = constants[40]*(exp((algebraic[110]+constants[135]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[127]+constants[139]*constants[95]*algebraic[71])/constants[140]))
    algebraic[195] = (algebraic[188]+algebraic[170]+algebraic[162])-algebraic[192]
    rates[47] = algebraic[195]
    algebraic[196] = algebraic[190]+algebraic[171]+algebraic[163]+algebraic[192]
    rates[48] = algebraic[196]
    algebraic[30] = custom_piecewise([greater(voi , constants[115]) & less(voi , constants[116]), constants[117] , True, constants[118]])
    algebraic[67] = constants[119]*algebraic[63]*constants[95]
    algebraic[194] = algebraic[185]*(constants[138]-constants[134])+algebraic[186]*(constants[139]-constants[135])+algebraic[188]*(constants[137]-constants[133])+algebraic[190]*(constants[138]-constants[134])+algebraic[192]*(constants[139]-constants[135])+algebraic[182]*(constants[137]-constants[133])+algebraic[171]*(constants[136]-constants[132])+algebraic[170]*(constants[136]-constants[132])+algebraic[169]*(constants[136]-constants[132])+algebraic[180]*(constants[139]-constants[135])+algebraic[178]*(constants[138]-constants[134])+algebraic[176]*(constants[137]-constants[133])+algebraic[174]*(constants[139]-constants[135])+algebraic[173]*(constants[138]-constants[134])+algebraic[172]*(constants[137]-constants[133])+algebraic[168]*(constants[136]-constants[132])+algebraic[167]*(constants[136]-constants[132])+algebraic[165]*(constants[136]-constants[132])
    algebraic[197] = constants[95]*((algebraic[194]-constants[119]*(algebraic[120]+algebraic[142]))-constants[120]*(algebraic[122]+algebraic[150]))
    algebraic[198] = (-algebraic[67]-algebraic[197])+algebraic[30]
    rates[0] = algebraic[198]
    return(rates)

def computeAlgebraic(constants, states, voi):
    algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
    states = array(states)
    voi = array(voi)
    algebraic[9] = constants[104]+states[10]+states[17]+states[21]+states[16]
    algebraic[31] = constants[121]*constants[122]*log(constants[66]*algebraic[9])
    algebraic[18] = constants[113]+states[22]
    algebraic[37] = constants[121]*constants[122]*log(constants[90]*algebraic[18])
    algebraic[22] = constants[114]+states[23]
    algebraic[43] = constants[121]*constants[122]*log(constants[91]*algebraic[22])
    algebraic[48] = constants[47]*(exp((algebraic[31]+algebraic[37])/(constants[121]*constants[122]))-exp(algebraic[43]/(constants[121]*constants[122])))
    algebraic[35] = states[0]/constants[94]
    algebraic[11] = constants[105]+states[11]+states[15]
    algebraic[42] = constants[121]*constants[122]*log(constants[67]*algebraic[11])
    algebraic[47] = constants[121]*constants[122]*log(constants[66]*algebraic[9])
    algebraic[60] = constants[121]*constants[122]*log(constants[74]*states[36])
    algebraic[63] = constants[14]*exp(algebraic[60]/(constants[121]*constants[122]))*(exp((algebraic[42]+constants[119]*constants[95]*algebraic[35])/(constants[121]*constants[122]))-exp(algebraic[47]/(constants[121]*constants[122])))
    algebraic[51] = constants[121]*constants[122]*log(constants[71]*states[33])
    algebraic[66] = constants[15]*(exp(algebraic[60]/(constants[121]*constants[122]))-exp((algebraic[51]+constants[130]*algebraic[47])/(constants[121]*constants[122])))
    algebraic[54] = constants[121]*constants[122]*log(constants[72]*states[34])
    algebraic[70] = constants[16]*(exp((algebraic[51]+constants[129]*algebraic[47])/(constants[121]*constants[122]))-exp(algebraic[54]/(constants[121]*constants[122])))
    algebraic[0] = constants[96]+states[1]
    algebraic[40] = constants[121]*constants[122]*log(constants[48]*algebraic[0])
    algebraic[8] = constants[97]+states[2]
    algebraic[45] = constants[121]*constants[122]*log(constants[49]*algebraic[8])
    algebraic[10] = constants[98]+states[3]
    algebraic[49] = constants[121]*constants[122]*log(constants[50]*algebraic[10])
    algebraic[68] = constants[0]*(exp((algebraic[40]+algebraic[45])/(constants[121]*constants[122]))-exp(algebraic[49]/(constants[121]*constants[122])))
    algebraic[20] = constants[101]+states[6]
    algebraic[58] = constants[121]*constants[122]*log(constants[53]*algebraic[20])
    algebraic[72] = constants[1]*(exp(algebraic[49]/(constants[121]*constants[122]))-exp((algebraic[58]+algebraic[45])/(constants[121]*constants[122])))
    algebraic[57] = constants[121]*constants[122]*log(constants[73]*states[35])
    algebraic[74] = constants[17]*(exp((algebraic[54]+constants[130]*algebraic[47])/(constants[121]*constants[122]))-exp(algebraic[57]/(constants[121]*constants[122])))
    algebraic[12] = constants[99]+states[4]
    algebraic[52] = constants[121]*constants[122]*log(constants[51]*algebraic[12])
    algebraic[16] = constants[100]+states[5]
    algebraic[55] = constants[121]*constants[122]*log(constants[52]*algebraic[16])
    algebraic[76] = constants[2]*(exp((algebraic[58]+algebraic[52])/(constants[121]*constants[122]))-exp(algebraic[55]/(constants[121]*constants[122])))
    algebraic[78] = constants[18]*(exp(algebraic[57]/(constants[121]*constants[122]))-exp((algebraic[60]+constants[129]*algebraic[47])/(constants[121]*constants[122])))
    algebraic[80] = constants[3]*(exp(algebraic[55]/(constants[121]*constants[122]))-exp((algebraic[40]+algebraic[52])/(constants[121]*constants[122])))
    algebraic[82] = ((constants[130]*algebraic[66]-constants[129]*algebraic[70])-constants[130]*algebraic[74])+constants[129]*algebraic[78]
    algebraic[25] = constants[102]+states[7]
    algebraic[61] = constants[121]*constants[122]*log(constants[54]*algebraic[25])
    algebraic[32] = constants[103]+states[8]
    algebraic[64] = constants[121]*constants[122]*log(constants[55]*algebraic[32])
    algebraic[84] = constants[4]*(exp((algebraic[52]+algebraic[61])/(constants[121]*constants[122]))-exp(algebraic[64]/(constants[121]*constants[122])))
    algebraic[73] = constants[121]*constants[122]*log(constants[64]*states[29])
    algebraic[87] = algebraic[73]
    algebraic[56] = constants[121]*constants[122]*log(constants[56]*states[24])
    algebraic[21] = constants[107]+states[14]
    algebraic[53] = constants[121]*constants[122]*log(constants[70]*algebraic[21])
    algebraic[89] = algebraic[56]+algebraic[53]
    algebraic[123] = constants[13]*(exp(algebraic[87]/(constants[121]*constants[122]))-exp(algebraic[89]/(constants[121]*constants[122])))
    algebraic[125] = algebraic[123]
    algebraic[13] = constants[109]+states[12]
    algebraic[50] = constants[121]*constants[122]*log(constants[68]*algebraic[13])
    algebraic[91] = algebraic[56]+algebraic[50]
    algebraic[59] = constants[121]*constants[122]*log(constants[57]*states[25])
    algebraic[93] = algebraic[59]
    algebraic[126] = constants[5]*(exp(algebraic[91]/(constants[121]*constants[122]))-exp(algebraic[93]/(constants[121]*constants[122])))
    algebraic[129] = -algebraic[126]
    algebraic[130] = algebraic[123]-algebraic[126]
    algebraic[69] = constants[121]*constants[122]*log(constants[60]*states[28])
    algebraic[95] = algebraic[69]
    algebraic[17] = constants[108]+states[13]
    algebraic[46] = constants[121]*constants[122]*log(constants[69]*algebraic[17])
    algebraic[77] = constants[121]*constants[122]*log(constants[61]*states[30])
    algebraic[97] = algebraic[46]+algebraic[77]
    algebraic[131] = constants[9]*(exp(algebraic[95]/(constants[121]*constants[122]))-exp(algebraic[97]/(constants[121]*constants[122])))
    algebraic[133] = algebraic[131]
    algebraic[1] = constants[106]+states[9]
    algebraic[41] = constants[121]*constants[122]*log(constants[65]*algebraic[1])
    algebraic[99] = algebraic[59]+algebraic[41]
    algebraic[62] = constants[121]*constants[122]*log(constants[58]*states[26])
    algebraic[101] = algebraic[62]
    algebraic[134] = constants[7]*(exp(algebraic[99]/(constants[121]*constants[122]))-exp(algebraic[101]/(constants[121]*constants[122])))
    algebraic[136] = algebraic[134]
    algebraic[29] = constants[121]*constants[122]*log(constants[66]*algebraic[9])
    algebraic[103] = algebraic[59]+constants[125]*algebraic[29]
    algebraic[65] = constants[121]*constants[122]*log(constants[59]*states[27])
    algebraic[105] = algebraic[65]
    algebraic[137] = constants[6]*(exp(algebraic[103]/(constants[121]*constants[122]))-exp(algebraic[105]/(constants[121]*constants[122])))
    algebraic[139] = -constants[125]*algebraic[137]
    algebraic[140] = algebraic[126]-algebraic[137]
    algebraic[107] = algebraic[65]
    algebraic[109] = algebraic[69]+constants[127]*algebraic[41]
    algebraic[141] = constants[8]*(exp(algebraic[107]/(constants[121]*constants[122]))-exp(algebraic[109]/(constants[121]*constants[122])))
    algebraic[143] = algebraic[137]-algebraic[141]
    algebraic[144] = algebraic[141]-algebraic[131]
    algebraic[111] = algebraic[77]
    algebraic[34] = constants[121]*constants[122]*log(constants[67]*algebraic[11])
    algebraic[81] = constants[121]*constants[122]*log(constants[62]*states[31])
    algebraic[113] = algebraic[81]+constants[126]*algebraic[34]
    algebraic[145] = constants[10]*(exp(algebraic[111]/(constants[121]*constants[122]))-exp(algebraic[113]/(constants[121]*constants[122])))
    algebraic[146] = constants[126]*algebraic[145]
    algebraic[147] = algebraic[131]-algebraic[145]
    algebraic[71] = states[0]/constants[94]
    algebraic[88] = constants[119]*constants[95]*algebraic[71]
    algebraic[75] = constants[140]*log(constants[66]*algebraic[9])
    algebraic[110] = constants[140]*log(constants[88]*states[47])
    algebraic[112] = algebraic[75]+algebraic[88]+algebraic[110]
    algebraic[15] = constants[110]+states[18]
    algebraic[79] = constants[140]*log(constants[75]*algebraic[15])
    algebraic[114] = algebraic[79]+algebraic[110]
    algebraic[120] = custom_piecewise([equal(algebraic[88] , 0.000000), constants[19]*(exp(algebraic[112]/constants[140])-exp(algebraic[114]/constants[140])) , True, (((constants[19]*algebraic[88])/constants[140])/(exp(algebraic[88]/constants[140])-1.00000))*(exp(algebraic[112]/constants[140])-exp(algebraic[114]/constants[140]))])
    algebraic[127] = constants[140]*log(constants[89]*states[48])
    algebraic[128] = algebraic[75]+algebraic[88]+algebraic[127]
    algebraic[132] = algebraic[79]+algebraic[127]
    algebraic[142] = custom_piecewise([equal(algebraic[88] , 0.000000), constants[20]*(exp(algebraic[128]/constants[140])-exp(algebraic[132]/constants[140])) , True, (((constants[20]*algebraic[88])/constants[140])/(exp(algebraic[88]/constants[140])-1.00000))*(exp(algebraic[128]/constants[140])-exp(algebraic[132]/constants[140]))])
    algebraic[149] = algebraic[142]+algebraic[120]
    algebraic[115] = algebraic[81]+constants[127]*algebraic[41]
    algebraic[85] = constants[121]*constants[122]*log(constants[63]*states[32])
    algebraic[117] = algebraic[85]
    algebraic[148] = constants[11]*(exp(algebraic[115]/(constants[121]*constants[122]))-exp(algebraic[117]/(constants[121]*constants[122])))
    algebraic[151] = algebraic[145]-algebraic[148]
    algebraic[90] = constants[120]*constants[95]*algebraic[71]
    algebraic[19] = constants[111]+states[19]
    algebraic[83] = constants[140]*log(constants[76]*algebraic[19])
    algebraic[116] = algebraic[83]+algebraic[90]+algebraic[110]
    algebraic[24] = constants[112]+states[20]
    algebraic[86] = constants[140]*log(constants[77]*algebraic[24])
    algebraic[118] = algebraic[86]+algebraic[110]
    algebraic[122] = custom_piecewise([equal(algebraic[90] , 0.000000), constants[21]*(exp(algebraic[116]/constants[140])-exp(algebraic[118]/constants[140])) , True, (((constants[21]*algebraic[90])/constants[140])/(exp(algebraic[90]/constants[140])-1.00000))*(exp(algebraic[116]/constants[140])-exp(algebraic[118]/constants[140]))])
    algebraic[135] = algebraic[83]+algebraic[90]+algebraic[127]
    algebraic[138] = algebraic[86]+algebraic[127]
    algebraic[150] = custom_piecewise([equal(algebraic[90] , 0.000000), constants[22]*(exp(algebraic[135]/constants[140])-exp(algebraic[138]/constants[140])) , True, (((constants[22]*algebraic[90])/constants[140])/(exp(algebraic[90]/constants[140])-1.00000))*(exp(algebraic[135]/constants[140])-exp(algebraic[138]/constants[140]))])
    algebraic[153] = -algebraic[122]-algebraic[150]
    algebraic[154] = algebraic[150]+algebraic[122]
    algebraic[119] = algebraic[85]
    algebraic[121] = algebraic[41]+algebraic[73]
    algebraic[152] = constants[12]*(exp(algebraic[119]/(constants[121]*constants[122]))-exp(algebraic[121]/(constants[121]*constants[122])))
    algebraic[156] = ((constants[127]*algebraic[141]-algebraic[134])-constants[127]*algebraic[148])+algebraic[152]
    algebraic[158] = algebraic[152]-algebraic[123]
    algebraic[157] = algebraic[148]-algebraic[152]
    algebraic[92] = constants[140]*log(constants[78]*states[37])
    algebraic[94] = constants[140]*log(constants[84]*states[43])
    algebraic[155] = constants[41]*(exp(algebraic[92]/constants[140])-exp((algebraic[94]+constants[131]*algebraic[75])/constants[140]))
    algebraic[96] = constants[140]*log(constants[79]*states[38])
    algebraic[98] = constants[140]*log(constants[85]*states[44])
    algebraic[159] = constants[43]*(exp(algebraic[96]/constants[140])-exp((algebraic[98]+constants[131]*algebraic[75])/constants[140]))
    algebraic[100] = constants[140]*log(constants[80]*states[39])
    algebraic[102] = constants[140]*log(constants[86]*states[45])
    algebraic[160] = constants[45]*(exp(algebraic[100]/constants[140])-exp((algebraic[102]+constants[131]*algebraic[75])/constants[140]))
    algebraic[104] = constants[140]*log(constants[81]*states[40])
    algebraic[106] = constants[140]*log(constants[87]*states[46])
    algebraic[161] = constants[42]*(exp(algebraic[104]/constants[140])-exp((algebraic[106]+constants[131]*algebraic[75])/constants[140]))
    algebraic[108] = constants[140]*log(constants[82]*states[41])
    algebraic[162] = constants[44]*(exp(algebraic[108]/constants[140])-exp((algebraic[110]+constants[131]*algebraic[75])/constants[140]))
    algebraic[124] = constants[140]*log(constants[83]*states[42])
    algebraic[163] = constants[46]*(exp(algebraic[124]/constants[140])-exp((algebraic[127]+constants[131]*algebraic[75])/constants[140]))
    algebraic[164] = constants[119]*(algebraic[161]+algebraic[162]+algebraic[160]+algebraic[159]+algebraic[155]+algebraic[163])
    algebraic[166] = (algebraic[164]-algebraic[120])-algebraic[142]
    algebraic[172] = constants[29]*(exp((algebraic[92]+constants[133]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[96]+constants[137]*constants[95]*algebraic[71])/constants[140]))
    algebraic[173] = constants[33]*(exp((algebraic[92]+constants[134]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[100]+constants[138]*constants[95]*algebraic[71])/constants[140]))
    algebraic[165] = constants[23]*(exp((algebraic[92]+constants[132]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[104]+constants[136]*constants[95]*algebraic[71])/constants[140]))
    algebraic[175] = ((-algebraic[165]-algebraic[172])-algebraic[173])-algebraic[155]
    algebraic[174] = constants[37]*(exp((algebraic[96]+constants[135]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[100]+constants[139]*constants[95]*algebraic[71])/constants[140]))
    algebraic[168] = constants[25]*(exp((algebraic[100]+constants[132]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[124]+constants[136]*constants[95]*algebraic[71])/constants[140]))
    algebraic[177] = (-algebraic[168]+algebraic[174]+algebraic[173])-algebraic[160]
    algebraic[176] = constants[30]*(exp((algebraic[104]+constants[133]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[108]+constants[137]*constants[95]*algebraic[71])/constants[140]))
    algebraic[167] = constants[24]*(exp((algebraic[96]+constants[132]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[108]+constants[136]*constants[95]*algebraic[71])/constants[140]))
    algebraic[179] = ((-algebraic[167]-algebraic[174])+algebraic[176])-algebraic[159]
    algebraic[178] = constants[34]*(exp((algebraic[104]+constants[134]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[124]+constants[138]*constants[95]*algebraic[71])/constants[140]))
    algebraic[181] = ((-algebraic[178]+algebraic[165])-algebraic[176])-algebraic[161]
    algebraic[180] = constants[38]*(exp((algebraic[108]+constants[135]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[124]+constants[139]*constants[95]*algebraic[71])/constants[140]))
    algebraic[183] = ((algebraic[176]+algebraic[167])-algebraic[162])-algebraic[180]
    algebraic[184] = ((algebraic[178]+algebraic[168])-algebraic[163])+algebraic[180]
    algebraic[169] = constants[26]*(exp((algebraic[94]+constants[132]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[106]+constants[136]*constants[95]*algebraic[71])/constants[140]))
    algebraic[182] = constants[31]*(exp((algebraic[94]+constants[133]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[98]+constants[137]*constants[95]*algebraic[71])/constants[140]))
    algebraic[185] = constants[35]*(exp((algebraic[94]+constants[134]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[102]+constants[138]*constants[95]*algebraic[71])/constants[140]))
    algebraic[187] = ((-algebraic[169]-algebraic[182])-algebraic[185])+algebraic[155]
    algebraic[171] = constants[28]*(exp((algebraic[102]+constants[132]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[127]+constants[136]*constants[95]*algebraic[71])/constants[140]))
    algebraic[186] = constants[39]*(exp((algebraic[98]+constants[135]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[102]+constants[139]*constants[95]*algebraic[71])/constants[140]))
    algebraic[189] = -algebraic[171]+algebraic[186]+algebraic[185]+algebraic[160]
    algebraic[170] = constants[27]*(exp((algebraic[98]+constants[132]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[110]+constants[136]*constants[95]*algebraic[71])/constants[140]))
    algebraic[188] = constants[32]*(exp((algebraic[106]+constants[133]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[110]+constants[137]*constants[95]*algebraic[71])/constants[140]))
    algebraic[191] = (-algebraic[170]-algebraic[186])+algebraic[188]+algebraic[159]
    algebraic[190] = constants[36]*(exp((algebraic[106]+constants[134]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[127]+constants[138]*constants[95]*algebraic[71])/constants[140]))
    algebraic[193] = ((-algebraic[190]+algebraic[169])-algebraic[188])+algebraic[161]
    algebraic[192] = constants[40]*(exp((algebraic[110]+constants[135]*constants[95]*algebraic[71])/constants[140])-exp((algebraic[127]+constants[139]*constants[95]*algebraic[71])/constants[140]))
    algebraic[195] = (algebraic[188]+algebraic[170]+algebraic[162])-algebraic[192]
    algebraic[196] = algebraic[190]+algebraic[171]+algebraic[163]+algebraic[192]
    algebraic[30] = custom_piecewise([greater(voi , constants[115]) & less(voi , constants[116]), constants[117] , True, constants[118]])
    algebraic[67] = constants[119]*algebraic[63]*constants[95]
    algebraic[194] = algebraic[185]*(constants[138]-constants[134])+algebraic[186]*(constants[139]-constants[135])+algebraic[188]*(constants[137]-constants[133])+algebraic[190]*(constants[138]-constants[134])+algebraic[192]*(constants[139]-constants[135])+algebraic[182]*(constants[137]-constants[133])+algebraic[171]*(constants[136]-constants[132])+algebraic[170]*(constants[136]-constants[132])+algebraic[169]*(constants[136]-constants[132])+algebraic[180]*(constants[139]-constants[135])+algebraic[178]*(constants[138]-constants[134])+algebraic[176]*(constants[137]-constants[133])+algebraic[174]*(constants[139]-constants[135])+algebraic[173]*(constants[138]-constants[134])+algebraic[172]*(constants[137]-constants[133])+algebraic[168]*(constants[136]-constants[132])+algebraic[167]*(constants[136]-constants[132])+algebraic[165]*(constants[136]-constants[132])
    algebraic[197] = constants[95]*((algebraic[194]-constants[119]*(algebraic[120]+algebraic[142]))-constants[120]*(algebraic[122]+algebraic[150]))
    algebraic[198] = (-algebraic[67]-algebraic[197])+algebraic[30]
    algebraic[2] = states[10]/constants[128]
    algebraic[3] = states[11]/constants[141]
    algebraic[4] = states[9]/constants[142]
    algebraic[5] = states[13]/constants[128]
    algebraic[6] = states[12]/constants[128]
    algebraic[7] = states[14]/constants[128]
    algebraic[14] = algebraic[8]+algebraic[10]
    algebraic[23] = algebraic[19]
    algebraic[26] = constants[125]*states[27]+constants[125]*states[28]+constants[126]*states[30]
    algebraic[27] = algebraic[0]+algebraic[10]+algebraic[20]+algebraic[16]
    algebraic[28] = algebraic[24]
    algebraic[33] = constants[130]*states[36]+constants[129]*states[34]+(constants[129]+constants[130])*states[35]
    algebraic[36] = algebraic[30]/(constants[119]*constants[95])
    algebraic[38] = constants[131]*(states[37]+states[38]+states[39]+states[40]+states[41]+states[42])
    algebraic[39] = algebraic[25]+algebraic[32]
    algebraic[44] = algebraic[9]+algebraic[22]+algebraic[11]+algebraic[15]+algebraic[38]+algebraic[26]+algebraic[33]
    return algebraic

def custom_piecewise(cases):
    """Compute result of a piecewise function"""
    return select(cases[0::2],cases[1::2])

def solve_model():
    (legend_states, legend_algebraic, legend_voi, legend_constants) = createLegends()
    """Solve model with ODE solver"""
    from scipy.integrate import ode
    # Initialise constants and state variables. Multiply vertain enzyme concentrations with Gmult factor
    glabels = [
        # 'q_Gs_init in component environment (fmol)',
    ]
    [igs, _, _] = find_indices(glabels, legend_constants, legend_states, legend_algebraic);
    Gmult = 1e0;  # 5e1;
    (init_states, constants) = initConsts(Gmult, igs)

    # Set timespan to solve over
    voi = linspace(0, 1e-7, 300)

    # options for stimulus protocol
    stimLabels = ["pulse_start in component environment (second)",
                  "pulse_end in component environment (second)",
                  "pulseMag in component environment (fA)",
                  "pulseHolding in component environment (fA)"
                  ]
    sd = find_indices_from_labels(stimLabels, legend_constants)
    constants[sd['pulse_start']] = 1.1e-8
    constants[sd['pulse_end']] = 9e-8
    constants[sd['pulseMag']] = 3e11
    constants[sd['pulseHolding']] = 0

    # Construct ODE object to solve
    r = ode(computeRates)
    r.set_integrator('vode', method='bdf', atol=1e-07, rtol=1e-07, max_step=1e-8)
    r.set_initial_value(init_states, voi[0])
    r.set_f_params(constants)

    # Solve model
    states = array([[0.0] * len(voi)] * sizeStates)
    states[:, 0] = init_states
    for (i, t) in enumerate(voi[1:]):
        if r.successful():
            r.integrate(t)
            states[:, i + 1] = r.y
        else:
            break

    # Compute algebraic variables
    algebraic = computeAlgebraic(constants, states, voi)
                                   

    for glab in glabels:
        print('Gmult applied for ' + glab)
    return (voi, constants, states, algebraic, Gmult)


def plot_model(voi, constants, states, algebraic, Gmult):
    """Plot variables against variable of integration"""

    (legend_states, legend_algebraic, legend_voi, legend_constants) = createLegends()

    labels = ['I_pulse in component environment (fA)',
              'V_m in component environment (volt)',
            'q_Cai in component environment (fmol)',
              'q_Ca_SR in component environment (fmol)',
              'q_Cao in component environment (fmol)',
              'q_Ca_TRPN in component environment (fmol)',
              'q_TRPN in component environment (fmol)']

    [_, i_st, i_alg] = find_indices(labels, legend_constants, legend_states, legend_algebraic)

    if i_st:
        plot_selected(i_st, voi, states, legend_voi, legend_states, ('QUANTITIES\nGmult=%g' % (Gmult)),
                  int(ceil(sqrt(len(i_st)))))
    if i_alg:
        plot_selected(i_alg, voi, algebraic, legend_voi, legend_algebraic, ('QUANTITIES\nGmult=%g' % (Gmult)),
                  int(ceil(sqrt(len(i_alg)))))

    # fluxes
    labels = ['v_CaGHK_i in component LCC (fmol_per_sec)',
              'v_CaGHK_o in component LCC (fmol_per_sec)',
              'v_RyR in component RyR (fmol_per_sec)',
              'v_Cai in component SERCA (fmol_per_sec)',
              'v_Ca_SR in component SERCA (fmol_per_sec)',
              'v_R_TRPNCa in component TRPN (fmol_per_sec)'
              ]
    [_, i_st, i_alg] = find_indices(labels, legend_constants, legend_states, legend_algebraic)
    if i_st:
        plot_selected(i_st, voi, states, legend_voi, legend_states, ('FLUXES\nGmult=%g' % (Gmult)),
                      int(ceil(sqrt(len(i_st)))))
    if i_alg:
        plot_selected(i_alg, voi, algebraic, legend_voi, legend_algebraic, ('FLUXES\nGmult=%g' % (Gmult)),
                      int(ceil(sqrt(len(i_alg)))))

    # plot all q variables to confirm if we are SS
    if False:
        plt.figure()
        n = ceil(sqrt(size(states, 2)))
        for i in range(1, size(states, 2)):
            plt.subplot(n, n, i + 1)
            plt.plot(voi, states[:, i])
            xlabel(legend_voi)
            title(legend_states[i])

    plt.show()


def find_indices(labels, legend_constants, legend_states, legend_algebraic):
    i_con = []
    i_st = []
    i_alg = []
    for lab in labels:
        try:
            i_con.append(legend_constants.index(lab))
        except:
            try:
                i_alg.append(legend_algebraic.index(lab))
            except:
                i_st.append(legend_states.index(lab))

    return (i_con, i_st, i_alg)


def find_indices_from_labels(labels, legend_list):
    # return indices for a name in a dict, with the key being the name of the label, andits value as the index
    sd = dict()
    i_found = [legend_list.index(lab) for lab in labels]
    sd = {labels[i].split(' ')[0]: i_found[i] for i in range(len(labels))}
    return sd


def plot_selected(ids, x, y, legend_x, legend_y, titlestr, ns):
    istart = 30;
    plt.figure();
    for i, id in enumerate(ids):
        plt.subplot(ns, ns, i + 1)
        plt.plot(x[istart:-1], y[id, istart:-1])
        plt.xlabel('time (s)')
        if 'FLUXES' in titlestr:
            plt.legend([legend_y[id].split(' (')[0]])
        else:
            plt.legend([legend_y[id].split(' ')[0]])
    # plt.tight_layout()
    plt.subplots_adjust(wspace=0.2, hspace=0.4)
    plt.suptitle(titlestr)


if __name__ == "__main__":
    t = time.time()
    (voi, constants,states, algebraic,Gmult) = solve_model()
    print('elapsed = ', round(time.time() - t,3), ' s')
    plot_model(voi, constants,states, algebraic, Gmult)