Generated Code

The following is matlab code generated by the CellML API from this CellML file. (Back to language selection)

The raw code is available.

function [VOI, STATES, ALGEBRAIC, CONSTANTS] = mainFunction()
    % This is the "main function".  In Matlab, things work best if you rename this function to match the filename.
   [VOI, STATES, ALGEBRAIC, CONSTANTS] = solveModel();
end

function [algebraicVariableCount] = getAlgebraicVariableCount() 
    % Used later when setting a global variable with the number of algebraic variables.
    % Note: This is not the "main method".  
    algebraicVariableCount =145;
end
% There are a total of 73 entries in each of the rate and state variable arrays.
% There are a total of 188 entries in the constant variable array.
%

function [VOI, STATES, ALGEBRAIC, CONSTANTS] = solveModel()
    % Create ALGEBRAIC of correct size
    global algebraicVariableCount;  algebraicVariableCount = getAlgebraicVariableCount();
    % Initialise constants and state variables
    [INIT_STATES, CONSTANTS] = initConsts;

    % Set timespan to solve over 
    tspan = [0, 10];

    % Set numerical accuracy options for ODE solver
    options = odeset('RelTol', 1e-06, 'AbsTol', 1e-06, 'MaxStep', 1);

    % Solve model with ODE solver
    [VOI, STATES] = ode15s(@(VOI, STATES)computeRates(VOI, STATES, CONSTANTS), tspan, INIT_STATES, options);

    % Compute algebraic variables
    [RATES, ALGEBRAIC] = computeRates(VOI, STATES, CONSTANTS);
    ALGEBRAIC = computeAlgebraic(ALGEBRAIC, CONSTANTS, STATES, VOI);

    % Plot state variables against variable of integration
    [LEGEND_STATES, LEGEND_ALGEBRAIC, LEGEND_VOI, LEGEND_CONSTANTS] = createLegends();
    figure();
    plot(VOI, STATES);
    xlabel(LEGEND_VOI);
    l = legend(LEGEND_STATES);
    set(l,'Interpreter','none');
end

function [LEGEND_STATES, LEGEND_ALGEBRAIC, LEGEND_VOI, LEGEND_CONSTANTS] = createLegends()
    LEGEND_STATES = ''; LEGEND_ALGEBRAIC = ''; LEGEND_VOI = ''; LEGEND_CONSTANTS = '';
    LEGEND_VOI = strpad('time in component Environment (second)');
    LEGEND_CONSTANTS(:,1) = strpad('F in component Parameters (sec_A_per_mol)');
    LEGEND_CONSTANTS(:,2) = strpad('vol_cyt in component Parameters (cubic_Meter)');
    LEGEND_CONSTANTS(:,3) = strpad('vol_md in component Parameters (cubic_Meter)');
    LEGEND_CONSTANTS(:,4) = strpad('vol_SR in component Parameters (cubic_Meter)');
    LEGEND_CONSTANTS(:,5) = strpad('vol_cyt_Ca in component Parameters (cubic_Meter)');
    LEGEND_CONSTANTS(:,6) = strpad('K_mr in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,7) = strpad('n_H in component Parameters (dimensionless)');
    LEGEND_CONSTANTS(:,8) = strpad('v_RyR in component Parameters (per_sec)');
    LEGEND_CONSTANTS(:,9) = strpad('K_a in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,10) = strpad('K_b in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,11) = strpad('K_c in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,12) = strpad('K_c_rate in component Parameters (per_sec)');
    LEGEND_CONSTANTS(:,13) = strpad('v_IP3R in component Parameters (per_sec)');
    LEGEND_CONSTANTS(:,14) = strpad('v_IP3 in component Parameters (per_sec)');
    LEGEND_CONSTANTS(:,15) = strpad('conc_IP3_ref in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,16) = strpad('alpha_4 in component Parameters (dimensionless)');
    LEGEND_CONSTANTS(:,17) = strpad('k4 in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,18) = strpad('I_r in component Parameters (per_sec)');
    LEGEND_CONSTANTS(:,19) = strpad('d1 in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,20) = strpad('a4 in component Parameters (per_M_sec)');
    LEGEND_CONSTANTS(:,21) = strpad('d4 in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,22) = strpad('a5 in component Parameters (per_M_sec)');
    LEGEND_CONSTANTS(:,23) = strpad('d5 in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,24) = strpad('G_CaSOC_max in component Parameters (siemens2)');
    LEGEND_CONSTANTS(:,25) = strpad('K_SOC in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,26) = strpad('delta_CT in component Parameters (volt)');
    LEGEND_CONSTANTS(:,27) = strpad('J_NaCl_max in component Parameters (mol_per_sec)');
    LEGEND_CONSTANTS(:,28) = strpad('J_KCl_max in component Parameters (mol_per_sec)');
    LEGEND_CONSTANTS(:,29) = strpad('J_NKCC_max in component Parameters (mol_per_sec)');
    LEGEND_CONSTANTS(:,30) = strpad('K_NKCC_Na in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,31) = strpad('K_NKCC_K in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,32) = strpad('K_NKCC_Cl1 in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,33) = strpad('K_NKCC_Cl2 in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,34) = strpad('k_on in component Parameters (per_M_sec)');
    LEGEND_CONSTANTS(:,35) = strpad('k_off in component Parameters (per_sec)');
    LEGEND_CONSTANTS(:,36) = strpad('alpha1 in component Parameters (per_sec)');
    LEGEND_CONSTANTS(:,37) = strpad('alpha2 in component Parameters (per_sec)');
    LEGEND_CONSTANTS(:,38) = strpad('alpha3 in component Parameters (per_sec)');
    LEGEND_CONSTANTS(:,39) = strpad('G_ClCa_max in component Parameters (siemens2)');
    LEGEND_CONSTANTS(:,40) = strpad('G_K_b in component Parameters (siemens2)');
    LEGEND_CONSTANTS(:,41) = strpad('G_Na_b in component Parameters (siemens2)');
    LEGEND_CONSTANTS(:,42) = strpad('G_Ca_b in component Parameters (siemens2)');
    LEGEND_CONSTANTS(:,43) = strpad('G_Cl_b in component Parameters (siemens2)');
    LEGEND_CONSTANTS(:,44) = strpad('G_Kir_max in component Parameters (siemens2)');
    LEGEND_CONSTANTS(:,45) = strpad('conc_K_ref in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,46) = strpad('G_Kv_max in component Parameters (siemens2)');
    LEGEND_CONSTANTS(:,47) = strpad('n_ATP in component Parameters (dimensionless)');
    LEGEND_CONSTANTS(:,48) = strpad('G_KATP_max in component Parameters (siemens2)');
    LEGEND_CONSTANTS(:,49) = strpad('G_KCa_max in component Parameters (siemens2)');
    LEGEND_CONSTANTS(:,50) = strpad('Tau_PF in component Parameters (second)');
    LEGEND_CONSTANTS(:,51) = strpad('Tau_PS in component Parameters (second)');
    LEGEND_CONSTANTS(:,52) = strpad('K_m_K in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,53) = strpad('K_m_Na_alpha1 in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,54) = strpad('K_m_Na_alpha2 in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,55) = strpad('I_NaCa_max in component Parameters (ampere)');
    LEGEND_CONSTANTS(:,56) = strpad('K_m_Cai in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,57) = strpad('K_m_Cao in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,58) = strpad('K_m_Nai in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,59) = strpad('K_m_Nao in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,60) = strpad('k_sat in component Parameters (dimensionless)');
    LEGEND_CONSTANTS(:,61) = strpad('gamma in component Parameters (dimensionless)');
    LEGEND_CONSTANTS(:,62) = strpad('K_mCa_act in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,63) = strpad('G_VONa_max in component Parameters (siemens2)');
    LEGEND_CONSTANTS(:,64) = strpad('Tau_m in component Parameters (second)');
    LEGEND_CONSTANTS(:,65) = strpad('Tau_h in component Parameters (second)');
    LEGEND_CONSTANTS(:,66) = strpad('K_m_CaP in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,67) = strpad('G_CaL_max in component Parameters (siemens2)');
    LEGEND_CONSTANTS(:,68) = strpad('conc_Bf_tot in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,69) = strpad('k_Bf_on in component Parameters (per_M_sec)');
    LEGEND_CONSTANTS(:,70) = strpad('k_Bf_off in component Parameters (per_sec)');
    LEGEND_CONSTANTS(:,71) = strpad('conc_Calseq_SR_tot in component Parameters (Molar)');
    LEGEND_CONSTANTS(:,72) = strpad('k_Calseq_on in component Parameters (per_M_sec)');
    LEGEND_CONSTANTS(:,73) = strpad('k_Calseq_off in component Parameters (per_sec)');
    LEGEND_CONSTANTS(:,74) = strpad('f_cyt in component Parameters (dimensionless)');
    LEGEND_CONSTANTS(:,75) = strpad('f_md in component Parameters (dimensionless)');
    LEGEND_CONSTANTS(:,76) = strpad('R in component Parameters (J_per_mol_K)');
    LEGEND_CONSTANTS(:,77) = strpad('T in component Parameters (kelvin)');
    LEGEND_CONSTANTS(:,78) = strpad('Area in component Parameters (sq_meter)');
    LEGEND_CONSTANTS(:,79) = strpad('L in component Parameters (meter)');
    LEGEND_CONSTANTS(:,80) = strpad('hcon in component Parameters (dimensionless)');
    LEGEND_CONSTANTS(:,81) = strpad('D_Ca in component Parameters (sq_m_per_sec)');
    LEGEND_CONSTANTS(:,177) = strpad('D_K in component Parameters (sq_m_per_sec)');
    LEGEND_CONSTANTS(:,178) = strpad('D_Na in component Parameters (sq_m_per_sec)');
    LEGEND_CONSTANTS(:,82) = strpad('D_Cl in component Parameters (sq_m_per_sec)');
    LEGEND_STATES(:,1) = strpad('conc_K_cyt in component Potassium_Conc_cyt (Molar)');
    LEGEND_ALGEBRAIC(:,38) = strpad('I_K_cyt_b in component Background_Currents (ampere)');
    LEGEND_ALGEBRAIC(:,88) = strpad('I_K_cyt_ir in component K_Inward_Rectifier_Current (ampere)');
    LEGEND_ALGEBRAIC(:,96) = strpad('I_K_cyt_ATP in component K_ATP_Current (ampere)');
    LEGEND_ALGEBRAIC(:,93) = strpad('I_K_cyt_v in component K_Delayed_Rectifier_Current (ampere)');
    LEGEND_ALGEBRAIC(:,99) = strpad('I_K_cyt_Ca in component K_Ca_Current (ampere)');
    LEGEND_ALGEBRAIC(:,141) = strpad('I_NaK_cyt_alpha1 in component NaK_alpha_Current (ampere)');
    LEGEND_ALGEBRAIC(:,78) = strpad('J_cyt_KCl in component Cl_flux_KCl (mol_per_sec)');
    LEGEND_ALGEBRAIC(:,80) = strpad('J_cyt_NKCC in component Cl_flux_NKCC (mol_per_sec)');
    LEGEND_ALGEBRAIC(:,116) = strpad('J_Kdiff in component Md_Cyt_Flux_K (mol_per_sec)');
    LEGEND_STATES(:,2) = strpad('conc_K_md in component Potassium_Conc_md (Molar)');
    LEGEND_ALGEBRAIC(:,42) = strpad('I_K_md_b in component Background_Currents (ampere)');
    LEGEND_ALGEBRAIC(:,91) = strpad('I_K_md_ir in component K_Inward_Rectifier_Current (ampere)');
    LEGEND_ALGEBRAIC(:,97) = strpad('I_K_md_ATP in component K_ATP_Current (ampere)');
    LEGEND_ALGEBRAIC(:,95) = strpad('I_K_md_v in component K_Delayed_Rectifier_Current (ampere)');
    LEGEND_ALGEBRAIC(:,101) = strpad('I_K_md_Ca in component K_Ca_Current (ampere)');
    LEGEND_ALGEBRAIC(:,142) = strpad('I_NaK_md_alpha2 in component NaK_alpha_Current (ampere)');
    LEGEND_ALGEBRAIC(:,79) = strpad('J_md_KCl in component Cl_flux_KCl (mol_per_sec)');
    LEGEND_ALGEBRAIC(:,81) = strpad('J_md_NKCC in component Cl_flux_NKCC (mol_per_sec)');
    LEGEND_STATES(:,3) = strpad('conc_Na_cyt in component Sodium_Conc_cyt (Molar)');
    LEGEND_ALGEBRAIC(:,44) = strpad('I_Na_cyt_b in component Background_Currents (ampere)');
    LEGEND_ALGEBRAIC(:,109) = strpad('I_VONa_cyt in component VONa_Current (ampere)');
    LEGEND_ALGEBRAIC(:,74) = strpad('I_Na_cyt_SOC in component SOC_Current_Na (ampere)');
    LEGEND_ALGEBRAIC(:,76) = strpad('J_cyt_NaCl in component Cl_flux_NaCl (mol_per_sec)');
    LEGEND_ALGEBRAIC(:,118) = strpad('J_Nadiff in component Md_Cyt_Flux_Na (mol_per_sec)');
    LEGEND_STATES(:,4) = strpad('conc_Na_md in component Sodium_Conc_md (Molar)');
    LEGEND_ALGEBRAIC(:,46) = strpad('I_Na_md_b in component Background_Currents (ampere)');
    LEGEND_ALGEBRAIC(:,110) = strpad('I_VONa_md in component VONa_Current (ampere)');
    LEGEND_ALGEBRAIC(:,75) = strpad('I_Na_md_SOC in component SOC_Current_Na (ampere)');
    LEGEND_ALGEBRAIC(:,77) = strpad('J_md_NaCl in component Cl_flux_NaCl (mol_per_sec)');
    LEGEND_ALGEBRAIC(:,108) = strpad('I_NaCa_md in component NCX_Current (ampere)');
    LEGEND_STATES(:,5) = strpad('conc_Cl_cyt in component Chloride_Conc_cyt (Molar)');
    LEGEND_ALGEBRAIC(:,48) = strpad('I_Cl_cyt_b in component Background_Currents (ampere)');
    LEGEND_ALGEBRAIC(:,82) = strpad('I_Cl_cyt_Ca in component ClCa_Current (ampere)');
    LEGEND_ALGEBRAIC(:,122) = strpad('J_Cldiff in component Md_Cyt_Flux_Cl (mol_per_sec)');
    LEGEND_STATES(:,6) = strpad('conc_Cl_md in component Chloride_Conc_md (Molar)');
    LEGEND_ALGEBRAIC(:,50) = strpad('I_Cl_md_b in component Background_Currents (ampere)');
    LEGEND_ALGEBRAIC(:,83) = strpad('I_Cl_md_Ca in component ClCa_Current (ampere)');
    LEGEND_STATES(:,7) = strpad('conc_Ca_cyt in component Calcium_Conc_cyt (Molar)');
    LEGEND_ALGEBRAIC(:,52) = strpad('I_Ca_cyt_b in component Background_Currents (ampere)');
    LEGEND_ALGEBRAIC(:,127) = strpad('I_Ca_cyt_P in component Ca_Pump_Current (ampere)');
    LEGEND_ALGEBRAIC(:,72) = strpad('I_Ca_cyt_SOC in component SOC_Current_Ca (ampere)');
    LEGEND_ALGEBRAIC(:,112) = strpad('I_Ca_cyt_L in component CaL_Current (ampere)');
    LEGEND_ALGEBRAIC(:,132) = strpad('I_cyt_SERCa_IP3R in component SERCA_Currents_IP3R (ampere)');
    LEGEND_ALGEBRAIC(:,56) = strpad('I_cyt_RyR in component RyR_Currents (ampere)');
    LEGEND_ALGEBRAIC(:,120) = strpad('J_Cadiff in component Md_Cyt_Flux_Ca (mol_per_sec)');
    LEGEND_ALGEBRAIC(:,64) = strpad('I_cyt_IP3R in component IP3R_Currents (ampere)');
    LEGEND_STATES(:,8) = strpad('conc_CaMN__ in component MLCK_Activation_CaM (Molar)');
    LEGEND_STATES(:,9) = strpad('conc_BfCa_cyt in component Bf_Kinetics (Molar)');
    LEGEND_ALGEBRAIC(:,134) = strpad('I_cyt_SERCa_RyR in component SERCA_Currents_RyR (ampere)');
    LEGEND_STATES(:,10) = strpad('conc_CaM_C_ in component MLCK_Activation_CaM (Molar)');
    LEGEND_STATES(:,11) = strpad('conc_CaM_CM in component MLCK_Activation_CaM (Molar)');
    LEGEND_STATES(:,12) = strpad('conc_CaMNCM in component MLCK_Activation_CaM (Molar)');
    LEGEND_STATES(:,13) = strpad('conc_CaM___ in component MLCK_Activation_CaM (Molar)');
    LEGEND_STATES(:,14) = strpad('conc_BfCa_md in component Bf_Kinetics (Molar)');
    LEGEND_STATES(:,15) = strpad('conc_Ca_md in component Calcium_Conc_md (Molar)');
    LEGEND_ALGEBRAIC(:,54) = strpad('I_Ca_md_b in component Background_Currents (ampere)');
    LEGEND_ALGEBRAIC(:,128) = strpad('I_Ca_md_P in component Ca_Pump_Current (ampere)');
    LEGEND_ALGEBRAIC(:,73) = strpad('I_Ca_md_SOC in component SOC_Current_Ca (ampere)');
    LEGEND_ALGEBRAIC(:,114) = strpad('I_Ca_md_L in component CaL_Current (ampere)');
    LEGEND_ALGEBRAIC(:,135) = strpad('I_md_SERCa_RyR in component SERCA_Currents_RyR (ampere)');
    LEGEND_ALGEBRAIC(:,58) = strpad('I_md_RyR in component RyR_Currents (ampere)');
    LEGEND_ALGEBRAIC(:,68) = strpad('I_md_IP3R in component IP3R_Currents (ampere)');
    LEGEND_ALGEBRAIC(:,133) = strpad('I_md_SERCa_IP3R in component SERCA_Currents_IP3R (ampere)');
    LEGEND_STATES(:,16) = strpad('conc_CaMN__ in component MLCK_Activation_CaM_md (Molar)');
    LEGEND_STATES(:,17) = strpad('conc_CaM_C_ in component MLCK_Activation_CaM_md (Molar)');
    LEGEND_STATES(:,18) = strpad('conc_CaM_CM in component MLCK_Activation_CaM_md (Molar)');
    LEGEND_STATES(:,19) = strpad('conc_CaMNCM in component MLCK_Activation_CaM_md (Molar)');
    LEGEND_STATES(:,20) = strpad('conc_Ca_IP3R in component Calcium_Conc_SR (Molar)');
    LEGEND_STATES(:,21) = strpad('conc_Ca_RyR in component Calcium_Conc_SR (Molar)');
    LEGEND_STATES(:,22) = strpad('conc_CalseqCa_IP3R in component Calseq_Kinetics (Molar)');
    LEGEND_STATES(:,23) = strpad('conc_CalseqCa_RyR in component Calseq_Kinetics (Molar)');
    LEGEND_STATES(:,24) = strpad('V_m_cyt in component Membrane_Voltage (volt)');
    LEGEND_STATES(:,25) = strpad('V_m_md in component Membrane_Voltage (volt)');
    LEGEND_ALGEBRAIC(:,35) = strpad('E_K_cyt in component Nernst_Potentials (volt)');
    LEGEND_ALGEBRAIC(:,40) = strpad('E_K_md in component Nernst_Potentials (volt)');
    LEGEND_ALGEBRAIC(:,43) = strpad('E_Na_cyt in component Nernst_Potentials (volt)');
    LEGEND_ALGEBRAIC(:,45) = strpad('E_Na_md in component Nernst_Potentials (volt)');
    LEGEND_ALGEBRAIC(:,47) = strpad('E_Cl_cyt in component Nernst_Potentials (volt)');
    LEGEND_ALGEBRAIC(:,49) = strpad('E_Cl_md in component Nernst_Potentials (volt)');
    LEGEND_ALGEBRAIC(:,51) = strpad('E_Ca_cyt in component Nernst_Potentials (volt)');
    LEGEND_ALGEBRAIC(:,53) = strpad('E_Ca_md in component Nernst_Potentials (volt)');
    LEGEND_CONSTANTS(:,83) = strpad('conc_K_out in component Nernst_Potentials (Molar)');
    LEGEND_CONSTANTS(:,84) = strpad('conc_Na_out in component Nernst_Potentials (Molar)');
    LEGEND_CONSTANTS(:,85) = strpad('conc_Cl_out in component Nernst_Potentials (Molar)');
    LEGEND_CONSTANTS(:,86) = strpad('conc_Ca_out in component Nernst_Potentials (Molar)');
    LEGEND_CONSTANTS(:,87) = strpad('z_K in component Nernst_Potentials (dimensionless)');
    LEGEND_CONSTANTS(:,88) = strpad('z_Na in component Nernst_Potentials (dimensionless)');
    LEGEND_CONSTANTS(:,89) = strpad('z_Cl in component Nernst_Potentials (dimensionless)');
    LEGEND_CONSTANTS(:,90) = strpad('z_Ca in component Nernst_Potentials (dimensionless)');
    LEGEND_CONSTANTS(:,179) = strpad('RTF in component Nernst_Potentials (volt)');
    LEGEND_ALGEBRAIC(:,131) = strpad('I_SERCA_max in component ROS_SERCA_Interaction (ampere)');
    LEGEND_ALGEBRAIC(:,125) = strpad('K_mf in component SERCA_Activation (Molar)');
    LEGEND_ALGEBRAIC(:,55) = strpad('P_cyt_RyR in component RyR_Currents (dimensionless)');
    LEGEND_ALGEBRAIC(:,57) = strpad('P_md_RyR in component RyR_Currents (dimensionless)');
    LEGEND_STATES(:,26) = strpad('w_cyt in component RyR_Currents (dimensionless)');
    LEGEND_STATES(:,27) = strpad('w_md in component RyR_Currents (dimensionless)');
    LEGEND_ALGEBRAIC(:,1) = strpad('w_inf_cyt in component RyR_Currents (dimensionless)');
    LEGEND_ALGEBRAIC(:,2) = strpad('w_inf_md in component RyR_Currents (dimensionless)');
    LEGEND_ALGEBRAIC(:,63) = strpad('x110_cyt in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_ALGEBRAIC(:,67) = strpad('x110_md in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_STATES(:,28) = strpad('conc_IP3_cyt in component IP3R_Currents (Molar)');
    LEGEND_STATES(:,29) = strpad('conc_IP3_md in component IP3R_Currents (Molar)');
    LEGEND_STATES(:,30) = strpad('x000_cyt in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_STATES(:,31) = strpad('x000_md in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_STATES(:,32) = strpad('x001_cyt in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_STATES(:,33) = strpad('x001_md in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_STATES(:,34) = strpad('x010_cyt in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_STATES(:,35) = strpad('x010_md in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_ALGEBRAIC(:,65) = strpad('x011_cyt in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_ALGEBRAIC(:,69) = strpad('x011_md in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_ALGEBRAIC(:,61) = strpad('x101_cyt in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_ALGEBRAIC(:,62) = strpad('x101_md in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_ALGEBRAIC(:,66) = strpad('x111_cyt in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_ALGEBRAIC(:,70) = strpad('x111_md in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_ALGEBRAIC(:,59) = strpad('x100_cyt in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_ALGEBRAIC(:,60) = strpad('x100_md in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_CONSTANTS(:,91) = strpad('d3 in component IP3R_Binding_Sites (Molar)');
    LEGEND_ALGEBRAIC(:,71) = strpad('conc_Ca_SR in component SOC_Current_Ca (Molar)');
    LEGEND_CONSTANTS(:,92) = strpad('z_Na in component SOC_Current_Na (dimensionless)');
    LEGEND_CONSTANTS(:,93) = strpad('z_Ca in component SOC_Current_Na (dimensionless)');
    LEGEND_CONSTANTS(:,94) = strpad('P_SOC_ratio in component SOC_Current_Na (dimensionless)');
    LEGEND_STATES(:,36) = strpad('y_C1_cyt in component Cl_Channels_cyt (dimensionless)');
    LEGEND_STATES(:,37) = strpad('y_C2_cyt in component Cl_Channels_cyt (dimensionless)');
    LEGEND_STATES(:,38) = strpad('y_C3_cyt in component Cl_Channels_cyt (dimensionless)');
    LEGEND_STATES(:,39) = strpad('y_O1_cyt in component Cl_Channels_cyt (dimensionless)');
    LEGEND_STATES(:,40) = strpad('y_O2_cyt in component Cl_Channels_cyt (dimensionless)');
    LEGEND_STATES(:,41) = strpad('y_O3_cyt in component Cl_Channels_cyt (dimensionless)');
    LEGEND_ALGEBRAIC(:,20) = strpad('Beta_cyt in component Cl_Channel_Rates (per_sec)');
    LEGEND_ALGEBRAIC(:,3) = strpad('y_C0_cyt in component Cl_Channels_cyt (dimensionless)');
    LEGEND_STATES(:,42) = strpad('y_C1_md in component Cl_Channels_md (dimensionless)');
    LEGEND_STATES(:,43) = strpad('y_C2_md in component Cl_Channels_md (dimensionless)');
    LEGEND_STATES(:,44) = strpad('y_C3_md in component Cl_Channels_md (dimensionless)');
    LEGEND_STATES(:,45) = strpad('y_O1_md in component Cl_Channels_md (dimensionless)');
    LEGEND_STATES(:,46) = strpad('y_O2_md in component Cl_Channels_md (dimensionless)');
    LEGEND_STATES(:,47) = strpad('y_O3_md in component Cl_Channels_md (dimensionless)');
    LEGEND_ALGEBRAIC(:,21) = strpad('Beta_md in component Cl_Channel_Rates (per_sec)');
    LEGEND_ALGEBRAIC(:,4) = strpad('y_C0_md in component Cl_Channels_md (dimensionless)');
    LEGEND_CONSTANTS(:,95) = strpad('V1 in component Cl_Channel_Rates (dimensionless)');
    LEGEND_CONSTANTS(:,96) = strpad('V2 in component Cl_Channel_Rates (per_V)');
    LEGEND_CONSTANTS(:,97) = strpad('Lambda_Beta in component Cl_Channel_Rates (per_sec)');
    LEGEND_ALGEBRAIC(:,87) = strpad('x_Kir_cyt in component K_Inward_Rectifier_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,90) = strpad('x_Kir_md in component K_Inward_Rectifier_Current (dimensionless)');
    LEGEND_ALGEBRAIC(:,84) = strpad('alpha_Kir_cyt in component K_ir_Rates (per_sec)');
    LEGEND_ALGEBRAIC(:,85) = strpad('alpha_Kir_md in component K_ir_Rates (per_sec)');
    LEGEND_ALGEBRAIC(:,86) = strpad('Beta_Kir_cyt in component K_ir_Rates (per_sec)');
    LEGEND_ALGEBRAIC(:,89) = strpad('Beta_Kir_md in component K_ir_Rates (per_sec)');
    LEGEND_ALGEBRAIC(:,92) = strpad('P_Kv_cyt in component Kv_Activations (dimensionless)');
    LEGEND_ALGEBRAIC(:,94) = strpad('P_Kv_md in component Kv_Activations (dimensionless)');
    LEGEND_STATES(:,48) = strpad('P1_cyt in component Kv_Activations (dimensionless)');
    LEGEND_STATES(:,49) = strpad('P1_md in component Kv_Activations (dimensionless)');
    LEGEND_STATES(:,50) = strpad('P2_cyt in component Kv_Activations (dimensionless)');
    LEGEND_STATES(:,51) = strpad('P2_md in component Kv_Activations (dimensionless)');
    LEGEND_ALGEBRAIC(:,5) = strpad('P_bar_Kv_cyt in component Kv_Activations (dimensionless)');
    LEGEND_ALGEBRAIC(:,6) = strpad('P_bar_Kv_md in component Kv_Activations (dimensionless)');
    LEGEND_ALGEBRAIC(:,22) = strpad('Tau_P1_cyt in component Kv_Time_Constants (second)');
    LEGEND_ALGEBRAIC(:,23) = strpad('Tau_P1_md in component Kv_Time_Constants (second)');
    LEGEND_ALGEBRAIC(:,24) = strpad('Tau_P2_cyt in component Kv_Time_Constants (second)');
    LEGEND_ALGEBRAIC(:,25) = strpad('Tau_P2_md in component Kv_Time_Constants (second)');
    LEGEND_ALGEBRAIC(:,98) = strpad('P_KCa_cyt in component KCa_Activations (dimensionless)');
    LEGEND_ALGEBRAIC(:,100) = strpad('P_KCa_md in component KCa_Activations (dimensionless)');
    LEGEND_STATES(:,52) = strpad('PF_cyt in component KCa_Activations (dimensionless)');
    LEGEND_STATES(:,53) = strpad('PF_md in component KCa_Activations (dimensionless)');
    LEGEND_STATES(:,54) = strpad('PS_cyt in component KCa_Activations (dimensionless)');
    LEGEND_STATES(:,55) = strpad('PS_md in component KCa_Activations (dimensionless)');
    LEGEND_ALGEBRAIC(:,26) = strpad('P_bar_KCa_cyt in component KCa_Activations (dimensionless)');
    LEGEND_ALGEBRAIC(:,27) = strpad('P_bar_KCa_md in component KCa_Activations (dimensionless)');
    LEGEND_ALGEBRAIC(:,7) = strpad('V_KCa_cyt in component cGMP_Ca_Interaction (volt)');
    LEGEND_ALGEBRAIC(:,8) = strpad('V_KCa_md in component cGMP_Ca_Interaction (volt)');
    LEGEND_ALGEBRAIC(:,139) = strpad('I_NaK_alpha1_max in component ROS_NaK_Interaction (ampere)');
    LEGEND_ALGEBRAIC(:,143) = strpad('I_NaK_alpha2_max in component NaK_alpha_Current (ampere)');
    LEGEND_ALGEBRAIC(:,102) = strpad('Psi_NaK_cyt in component NaK_alpha_Waveform (dimensionless)');
    LEGEND_ALGEBRAIC(:,103) = strpad('Psi_NaK_md in component NaK_alpha_Waveform (dimensionless)');
    LEGEND_CONSTANTS(:,180) = strpad('Sigma in component NaK_alpha_Waveform (dimensionless)');
    LEGEND_ALGEBRAIC(:,104) = strpad('K_a_NaCa in component NCX_Waveform (dimensionless)');
    LEGEND_ALGEBRAIC(:,105) = strpad('Psi_F in component NCX_Waveform (dimensionless)');
    LEGEND_ALGEBRAIC(:,106) = strpad('Psi_R in component NCX_Waveform (dimensionless)');
    LEGEND_ALGEBRAIC(:,107) = strpad('G in component NCX_Waveform (M_fourpow)');
    LEGEND_STATES(:,56) = strpad('m_VONa_cyt in component VONa_Channels (dimensionless)');
    LEGEND_STATES(:,57) = strpad('m_VONa_md in component VONa_Channels (dimensionless)');
    LEGEND_STATES(:,58) = strpad('h_VONa_cyt in component VONa_Channels (dimensionless)');
    LEGEND_STATES(:,59) = strpad('h_VONa_md in component VONa_Channels (dimensionless)');
    LEGEND_ALGEBRAIC(:,9) = strpad('m_bar_cyt in component VONa_Channels (dimensionless)');
    LEGEND_ALGEBRAIC(:,10) = strpad('m_bar_md in component VONa_Channels (dimensionless)');
    LEGEND_ALGEBRAIC(:,11) = strpad('h_bar_cyt in component VONa_Channels (dimensionless)');
    LEGEND_ALGEBRAIC(:,12) = strpad('h_bar_md in component VONa_Channels (dimensionless)');
    LEGEND_ALGEBRAIC(:,126) = strpad('I_CaP_max in component CaP_Current_Max (ampere)');
    LEGEND_STATES(:,60) = strpad('d_L_cyt in component CaL_Activations (dimensionless)');
    LEGEND_STATES(:,61) = strpad('d_L_md in component CaL_Activations (dimensionless)');
    LEGEND_ALGEBRAIC(:,111) = strpad('f_L_cyt in component CaL_Activations (dimensionless)');
    LEGEND_ALGEBRAIC(:,113) = strpad('f_L_md in component CaL_Activations (dimensionless)');
    LEGEND_STATES(:,62) = strpad('f_F_cyt in component CaL_Activations (dimensionless)');
    LEGEND_STATES(:,63) = strpad('f_F_md in component CaL_Activations (dimensionless)');
    LEGEND_ALGEBRAIC(:,13) = strpad('d_bar_L_cyt in component CaL_Activations (dimensionless)');
    LEGEND_ALGEBRAIC(:,14) = strpad('d_bar_L_md in component CaL_Activations (dimensionless)');
    LEGEND_ALGEBRAIC(:,15) = strpad('f_bar_F_cyt in component CaL_Activations (dimensionless)');
    LEGEND_ALGEBRAIC(:,16) = strpad('f_bar_F_md in component CaL_Activations (dimensionless)');
    LEGEND_ALGEBRAIC(:,28) = strpad('Tau_d_cyt in component CaL_Time_Constants (second)');
    LEGEND_ALGEBRAIC(:,29) = strpad('Tau_d_md in component CaL_Time_Constants (second)');
    LEGEND_ALGEBRAIC(:,30) = strpad('Tau_f_cyt in component CaL_Time_Constants (second)');
    LEGEND_ALGEBRAIC(:,31) = strpad('Tau_f_md in component CaL_Time_Constants (second)');
    LEGEND_ALGEBRAIC(:,115) = strpad('Xi_K in component Md_Cyt_Flux_K (dimensionless)');
    LEGEND_CONSTANTS(:,183) = strpad('P_K_diff in component Md_Cyt_Flux_K (m_per_sec)');
    LEGEND_CONSTANTS(:,98) = strpad('z_K in component Md_Cyt_Flux_K (dimensionless)');
    LEGEND_CONSTANTS(:,187) = strpad('Volsec in component Md_Cyt_Flux_K (cubic_m_per_sec)');
    LEGEND_CONSTANTS(:,99) = strpad('DiffArea in component Md_Cyt_Flux_K (sq_meter)');
    LEGEND_ALGEBRAIC(:,117) = strpad('Xi_Na in component Md_Cyt_Flux_Na (dimensionless)');
    LEGEND_CONSTANTS(:,184) = strpad('P_Na_diff in component Md_Cyt_Flux_Na (m_per_sec)');
    LEGEND_CONSTANTS(:,100) = strpad('z_Na in component Md_Cyt_Flux_Na (dimensionless)');
    LEGEND_CONSTANTS(:,188) = strpad('Volsec in component Md_Cyt_Flux_Na (cubic_m_per_sec)');
    LEGEND_CONSTANTS(:,101) = strpad('DiffArea in component Md_Cyt_Flux_Na (sq_meter)');
    LEGEND_ALGEBRAIC(:,119) = strpad('Xi_Ca in component Md_Cyt_Flux_Ca (dimensionless)');
    LEGEND_CONSTANTS(:,181) = strpad('P_Ca_diff in component Md_Cyt_Flux_Ca (m_per_sec)');
    LEGEND_CONSTANTS(:,102) = strpad('z_Ca in component Md_Cyt_Flux_Ca (dimensionless)');
    LEGEND_CONSTANTS(:,185) = strpad('Volsec in component Md_Cyt_Flux_Ca (cubic_m_per_sec)');
    LEGEND_CONSTANTS(:,103) = strpad('DiffArea in component Md_Cyt_Flux_Ca (sq_meter)');
    LEGEND_ALGEBRAIC(:,121) = strpad('Xi_Cl in component Md_Cyt_Flux_Cl (dimensionless)');
    LEGEND_CONSTANTS(:,182) = strpad('P_Cl_diff in component Md_Cyt_Flux_Cl (m_per_sec)');
    LEGEND_CONSTANTS(:,104) = strpad('z_Cl in component Md_Cyt_Flux_Cl (dimensionless)');
    LEGEND_CONSTANTS(:,186) = strpad('Volsec in component Md_Cyt_Flux_Cl (cubic_m_per_sec)');
    LEGEND_CONSTANTS(:,105) = strpad('DiffArea in component Md_Cyt_Flux_Cl (sq_meter)');
    LEGEND_ALGEBRAIC(:,144) = strpad('I_cyt_tot in component Membrane_Potential_Cyt (ampere)');
    LEGEND_ALGEBRAIC(:,136) = strpad('I_cyt_up in component Membrane_Potential_Cyt (ampere)');
    LEGEND_ALGEBRAIC(:,145) = strpad('I_md_tot in component Membrane_Potential_md (ampere)');
    LEGEND_CONSTANTS(:,106) = strpad('C_m in component Membrane_Voltage (farad2)');
    LEGEND_ALGEBRAIC(:,138) = strpad('I_md_up in component Membrane_Potential_md (ampere)');
    LEGEND_CONSTANTS(:,107) = strpad('k_NO_O2 in component Parameters_VR (per_M_sec)');
    LEGEND_CONSTANTS(:,108) = strpad('k_NO_cons in component Parameters_VR (per_sec)');
    LEGEND_CONSTANTS(:,109) = strpad('k_SOD_O2 in component Parameters_VR (per_M_sec)');
    LEGEND_CONSTANTS(:,110) = strpad('k_cat in component Parameters_VR (per_M_sec)');
    LEGEND_CONSTANTS(:,111) = strpad('conc_SOD in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,112) = strpad('conc_CAT in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,113) = strpad('J_NO in component Parameters_VR (mol_per_sq_m_sec)');
    LEGEND_CONSTANTS(:,114) = strpad('G_O2 in component Parameters_VR (M_per_sec)');
    LEGEND_CONSTANTS(:,115) = strpad('J_H2O2 in component Parameters_VR (mol_per_sq_m_sec)');
    LEGEND_CONSTANTS(:,116) = strpad('P_NO in component Parameters_VR (m_per_sec)');
    LEGEND_CONSTANTS(:,117) = strpad('P_O2 in component Parameters_VR (m_per_sec)');
    LEGEND_CONSTANTS(:,118) = strpad('P_H2O2 in component Parameters_VR (m_per_sec)');
    LEGEND_CONSTANTS(:,119) = strpad('k_sGC_for in component Parameters_VR (per_M_sec)');
    LEGEND_CONSTANTS(:,120) = strpad('k_sGC_back in component Parameters_VR (per_sec)');
    LEGEND_CONSTANTS(:,121) = strpad('k2_sGC_for in component Parameters_VR (per_sec)');
    LEGEND_CONSTANTS(:,122) = strpad('k3_sCG_for in component Parameters_VR (per_M_sec)');
    LEGEND_CONSTANTS(:,123) = strpad('K4_sGC_for in component Parameters_VR (per_M_sec)');
    LEGEND_CONSTANTS(:,124) = strpad('V_sGC_max in component Parameters_VR (M_per_sec)');
    LEGEND_CONSTANTS(:,125) = strpad('k_PDE in component Parameters_VR (per_sec)');
    LEGEND_CONSTANTS(:,126) = strpad('K_M_PDE in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,127) = strpad('K_MLCP in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,128) = strpad('K_SERCA in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,129) = strpad('K_PMCA in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,130) = strpad('K_KCa_cGMP in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,131) = strpad('V_KCa_cGMP in component Parameters_VR (volt)');
    LEGEND_CONSTANTS(:,132) = strpad('K_KCa_NO in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,133) = strpad('V_KCa_NO in component Parameters_VR (volt)');
    LEGEND_CONSTANTS(:,134) = strpad('V_Ca in component Parameters_VR (volt)');
    LEGEND_CONSTANTS(:,135) = strpad('V_B in component Parameters_VR (volt)');
    LEGEND_CONSTANTS(:,136) = strpad('k_KCa in component Parameters_VR (volt)');
    LEGEND_CONSTANTS(:,137) = strpad('K_SERCA_O2 in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,138) = strpad('K_SERCA_H2O2 in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,139) = strpad('K_NaK_O2 in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,140) = strpad('K_NaK_H202 in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,141) = strpad('K_MLCP_H2O2 in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,142) = strpad('alpha_MLCP in component Parameters_VR (dimensionless)');
    LEGEND_CONSTANTS(:,143) = strpad('k1_CAM_on in component Parameters_VR (per_sq_M_sec)');
    LEGEND_CONSTANTS(:,144) = strpad('k1_CAM_off in component Parameters_VR (per_sec)');
    LEGEND_CONSTANTS(:,145) = strpad('k2_CAM_on in component Parameters_VR (per_sq_M_sec)');
    LEGEND_CONSTANTS(:,146) = strpad('k2_CAM_off in component Parameters_VR (per_sec)');
    LEGEND_CONSTANTS(:,147) = strpad('k3_CAM_on in component Parameters_VR (per_sq_M_sec)');
    LEGEND_CONSTANTS(:,148) = strpad('k3_CAM_off in component Parameters_VR (per_sec)');
    LEGEND_CONSTANTS(:,149) = strpad('k4_CAM_on in component Parameters_VR (per_sq_M_sec)');
    LEGEND_CONSTANTS(:,150) = strpad('k4_CAM_off in component Parameters_VR (per_sec)');
    LEGEND_CONSTANTS(:,151) = strpad('k5_CAM_on in component Parameters_VR (per_M_sec)');
    LEGEND_CONSTANTS(:,152) = strpad('k5_CAM_off in component Parameters_VR (per_sec)');
    LEGEND_CONSTANTS(:,153) = strpad('k6_CAM_on in component Parameters_VR (per_sq_M_sec)');
    LEGEND_CONSTANTS(:,154) = strpad('k6_CAM_off in component Parameters_VR (per_sec)');
    LEGEND_CONSTANTS(:,155) = strpad('k7_CAM_on in component Parameters_VR (per_M_sec)');
    LEGEND_CONSTANTS(:,156) = strpad('k7_CAM_off in component Parameters_VR (per_sec)');
    LEGEND_CONSTANTS(:,157) = strpad('k_Myo_MLCK in component Parameters_VR (per_M_sec)');
    LEGEND_CONSTANTS(:,158) = strpad('k_Myo_MLCP in component Parameters_VR (per_sec)');
    LEGEND_CONSTANTS(:,159) = strpad('k3_Myo in component Parameters_VR (per_sec)');
    LEGEND_CONSTANTS(:,160) = strpad('k4_Myo in component Parameters_VR (per_sec)');
    LEGEND_CONSTANTS(:,161) = strpad('k7_Myo in component Parameters_VR (per_sec)');
    LEGEND_CONSTANTS(:,162) = strpad('conc_CaM_tot in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,163) = strpad('conc_MLCK_tot in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,164) = strpad('conc_Myo_tot in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,165) = strpad('K_mf_rest in component Parameters_VR (Molar)');
    LEGEND_CONSTANTS(:,166) = strpad('I_PMCA_rest_max in component Parameters_VR (ampere)');
    LEGEND_CONSTANTS(:,167) = strpad('I_SERCA_rest_max in component Parameters_VR (ampere)');
    LEGEND_CONSTANTS(:,168) = strpad('I_NaK_alpha1_rest_max in component Parameters_VR (ampere)');
    LEGEND_CONSTANTS(:,169) = strpad('I_NaK_alpha2_rest_max in component Parameters_VR (ampere)');
    LEGEND_ALGEBRAIC(:,123) = strpad('conc_CaMNC_ in component MLCK_Activation_CaM (Molar)');
    LEGEND_ALGEBRAIC(:,124) = strpad('conc_MLCK_free in component MLCK_Activation_CaM (Molar)');
    LEGEND_STATES(:,64) = strpad('conc_Myo in component MLCK_Phospho_Myosin (Molar)');
    LEGEND_STATES(:,65) = strpad('conc_MyoP in component MLCK_Phospho_Myosin (Molar)');
    LEGEND_STATES(:,66) = strpad('conc_AMyo in component MLCK_Phospho_Myosin (Molar)');
    LEGEND_ALGEBRAIC(:,17) = strpad('conc_AMyoP in component MLCK_Phospho_Myosin (Molar)');
    LEGEND_ALGEBRAIC(:,32) = strpad('k1_6_Myo in component MLCP_Rate_Constants (per_sec)');
    LEGEND_ALGEBRAIC(:,39) = strpad('k2_Myo in component MLCP_Rate_Constants (per_sec)');
    LEGEND_ALGEBRAIC(:,41) = strpad('k5_Myo in component MLCP_Rate_Constants (per_sec)');
    LEGEND_CONSTANTS(:,170) = strpad('K_M_cGMP in component MLCP_Rate_Constants (Molar)');
    LEGEND_STATES(:,67) = strpad('conc_cGMP in component cGMP_Kinetics (Molar)');
    LEGEND_ALGEBRAIC(:,36) = strpad('R_MLCP in component MLCP_Rate_Constants (dimensionless)');
    LEGEND_STATES(:,68) = strpad('conc_H2O2_cyt in component H2O2_Kinetics (Molar)');
    LEGEND_STATES(:,69) = strpad('conc_NO_cyt in component NO_Kinetics (Molar)');
    LEGEND_STATES(:,70) = strpad('conc_O2_cyt in component O2_Kinetics (Molar)');
    LEGEND_ALGEBRAIC(:,18) = strpad('zeta in component Zeta (dimensionless)');
    LEGEND_STATES(:,71) = strpad('conc_Eb in component sGC_Activation (dimensionless)');
    LEGEND_ALGEBRAIC(:,19) = strpad('conc_E6c in component sGC_Activation (dimensionless)');
    LEGEND_STATES(:,72) = strpad('conc_E5c in component sGC_Activation (dimensionless)');
    LEGEND_ALGEBRAIC(:,33) = strpad('k4_sGC_for in component sGC_Activation (per_sec)');
    LEGEND_CONSTANTS(:,171) = strpad('convert in component cGMP_Ca_Interaction (dimensionless)');
    LEGEND_CONSTANTS(:,172) = strpad('K_mod in component SERCA_Activation (dimensionless)');
    LEGEND_ALGEBRAIC(:,129) = strpad('theta_SERCA_O2 in component ROS_SERCA_Interaction (dimensionless)');
    LEGEND_ALGEBRAIC(:,130) = strpad('theta_SERCA_H2O2 in component ROS_SERCA_Interaction (dimensionless)');
    LEGEND_STATES(:,73) = strpad('conc_CaM___ in component MLCK_Activation_CaM_md (Molar)');
    LEGEND_ALGEBRAIC(:,137) = strpad('conc_CaMNC_ in component MLCK_Activation_CaM_md (Molar)');
    LEGEND_ALGEBRAIC(:,140) = strpad('conc_MLCK_free in component MLCK_Activation_CaM_md (Molar)');
    LEGEND_ALGEBRAIC(:,34) = strpad('f_p in component Contractile_Force (dimensionless)');
    LEGEND_ALGEBRAIC(:,37) = strpad('F_contract in component Contractile_Force (dimensionless)');
    LEGEND_CONSTANTS(:,173) = strpad('a in component Contractile_Force (dimensionless)');
    LEGEND_CONSTANTS(:,174) = strpad('b in component Contractile_Force (dimensionless)');
    LEGEND_CONSTANTS(:,175) = strpad('c in component Contractile_Force (dimensionless)');
    LEGEND_CONSTANTS(:,176) = strpad('d in component Contractile_Force (dimensionless)');
    LEGEND_RATES(:,1) = strpad('d/dt conc_K_cyt in component Potassium_Conc_cyt (Molar)');
    LEGEND_RATES(:,2) = strpad('d/dt conc_K_md in component Potassium_Conc_md (Molar)');
    LEGEND_RATES(:,3) = strpad('d/dt conc_Na_cyt in component Sodium_Conc_cyt (Molar)');
    LEGEND_RATES(:,4) = strpad('d/dt conc_Na_md in component Sodium_Conc_md (Molar)');
    LEGEND_RATES(:,5) = strpad('d/dt conc_Cl_cyt in component Chloride_Conc_cyt (Molar)');
    LEGEND_RATES(:,6) = strpad('d/dt conc_Cl_md in component Chloride_Conc_md (Molar)');
    LEGEND_RATES(:,7) = strpad('d/dt conc_Ca_cyt in component Calcium_Conc_cyt (Molar)');
    LEGEND_RATES(:,8) = strpad('d/dt conc_CaMN__ in component MLCK_Activation_CaM (Molar)');
    LEGEND_RATES(:,10) = strpad('d/dt conc_CaM_C_ in component MLCK_Activation_CaM (Molar)');
    LEGEND_RATES(:,11) = strpad('d/dt conc_CaM_CM in component MLCK_Activation_CaM (Molar)');
    LEGEND_RATES(:,9) = strpad('d/dt conc_BfCa_cyt in component Bf_Kinetics (Molar)');
    LEGEND_RATES(:,14) = strpad('d/dt conc_BfCa_md in component Bf_Kinetics (Molar)');
    LEGEND_RATES(:,15) = strpad('d/dt conc_Ca_md in component Calcium_Conc_md (Molar)');
    LEGEND_RATES(:,16) = strpad('d/dt conc_CaMN__ in component MLCK_Activation_CaM_md (Molar)');
    LEGEND_RATES(:,17) = strpad('d/dt conc_CaM_C_ in component MLCK_Activation_CaM_md (Molar)');
    LEGEND_RATES(:,18) = strpad('d/dt conc_CaM_CM in component MLCK_Activation_CaM_md (Molar)');
    LEGEND_RATES(:,20) = strpad('d/dt conc_Ca_IP3R in component Calcium_Conc_SR (Molar)');
    LEGEND_RATES(:,22) = strpad('d/dt conc_CalseqCa_IP3R in component Calseq_Kinetics (Molar)');
    LEGEND_RATES(:,21) = strpad('d/dt conc_Ca_RyR in component Calcium_Conc_SR (Molar)');
    LEGEND_RATES(:,23) = strpad('d/dt conc_CalseqCa_RyR in component Calseq_Kinetics (Molar)');
    LEGEND_RATES(:,26) = strpad('d/dt w_cyt in component RyR_Currents (dimensionless)');
    LEGEND_RATES(:,27) = strpad('d/dt w_md in component RyR_Currents (dimensionless)');
    LEGEND_RATES(:,28) = strpad('d/dt conc_IP3_cyt in component IP3R_Currents (Molar)');
    LEGEND_RATES(:,29) = strpad('d/dt conc_IP3_md in component IP3R_Currents (Molar)');
    LEGEND_RATES(:,30) = strpad('d/dt x000_cyt in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_RATES(:,31) = strpad('d/dt x000_md in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_RATES(:,32) = strpad('d/dt x001_cyt in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_RATES(:,33) = strpad('d/dt x001_md in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_RATES(:,34) = strpad('d/dt x010_cyt in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_RATES(:,35) = strpad('d/dt x010_md in component IP3R_Binding_Sites (dimensionless)');
    LEGEND_RATES(:,36) = strpad('d/dt y_C1_cyt in component Cl_Channels_cyt (dimensionless)');
    LEGEND_RATES(:,37) = strpad('d/dt y_C2_cyt in component Cl_Channels_cyt (dimensionless)');
    LEGEND_RATES(:,38) = strpad('d/dt y_C3_cyt in component Cl_Channels_cyt (dimensionless)');
    LEGEND_RATES(:,39) = strpad('d/dt y_O1_cyt in component Cl_Channels_cyt (dimensionless)');
    LEGEND_RATES(:,40) = strpad('d/dt y_O2_cyt in component Cl_Channels_cyt (dimensionless)');
    LEGEND_RATES(:,41) = strpad('d/dt y_O3_cyt in component Cl_Channels_cyt (dimensionless)');
    LEGEND_RATES(:,42) = strpad('d/dt y_C1_md in component Cl_Channels_md (dimensionless)');
    LEGEND_RATES(:,43) = strpad('d/dt y_C2_md in component Cl_Channels_md (dimensionless)');
    LEGEND_RATES(:,44) = strpad('d/dt y_C3_md in component Cl_Channels_md (dimensionless)');
    LEGEND_RATES(:,45) = strpad('d/dt y_O1_md in component Cl_Channels_md (dimensionless)');
    LEGEND_RATES(:,46) = strpad('d/dt y_O2_md in component Cl_Channels_md (dimensionless)');
    LEGEND_RATES(:,47) = strpad('d/dt y_O3_md in component Cl_Channels_md (dimensionless)');
    LEGEND_RATES(:,48) = strpad('d/dt P1_cyt in component Kv_Activations (dimensionless)');
    LEGEND_RATES(:,49) = strpad('d/dt P1_md in component Kv_Activations (dimensionless)');
    LEGEND_RATES(:,50) = strpad('d/dt P2_cyt in component Kv_Activations (dimensionless)');
    LEGEND_RATES(:,51) = strpad('d/dt P2_md in component Kv_Activations (dimensionless)');
    LEGEND_RATES(:,52) = strpad('d/dt PF_cyt in component KCa_Activations (dimensionless)');
    LEGEND_RATES(:,53) = strpad('d/dt PF_md in component KCa_Activations (dimensionless)');
    LEGEND_RATES(:,54) = strpad('d/dt PS_cyt in component KCa_Activations (dimensionless)');
    LEGEND_RATES(:,55) = strpad('d/dt PS_md in component KCa_Activations (dimensionless)');
    LEGEND_RATES(:,56) = strpad('d/dt m_VONa_cyt in component VONa_Channels (dimensionless)');
    LEGEND_RATES(:,57) = strpad('d/dt m_VONa_md in component VONa_Channels (dimensionless)');
    LEGEND_RATES(:,58) = strpad('d/dt h_VONa_cyt in component VONa_Channels (dimensionless)');
    LEGEND_RATES(:,59) = strpad('d/dt h_VONa_md in component VONa_Channels (dimensionless)');
    LEGEND_RATES(:,60) = strpad('d/dt d_L_cyt in component CaL_Activations (dimensionless)');
    LEGEND_RATES(:,61) = strpad('d/dt d_L_md in component CaL_Activations (dimensionless)');
    LEGEND_RATES(:,62) = strpad('d/dt f_F_cyt in component CaL_Activations (dimensionless)');
    LEGEND_RATES(:,63) = strpad('d/dt f_F_md in component CaL_Activations (dimensionless)');
    LEGEND_RATES(:,24) = strpad('d/dt V_m_cyt in component Membrane_Voltage (volt)');
    LEGEND_RATES(:,25) = strpad('d/dt V_m_md in component Membrane_Voltage (volt)');
    LEGEND_RATES(:,13) = strpad('d/dt conc_CaM___ in component MLCK_Activation_CaM (Molar)');
    LEGEND_RATES(:,12) = strpad('d/dt conc_CaMNCM in component MLCK_Activation_CaM (Molar)');
    LEGEND_RATES(:,64) = strpad('d/dt conc_Myo in component MLCK_Phospho_Myosin (Molar)');
    LEGEND_RATES(:,65) = strpad('d/dt conc_MyoP in component MLCK_Phospho_Myosin (Molar)');
    LEGEND_RATES(:,66) = strpad('d/dt conc_AMyo in component MLCK_Phospho_Myosin (Molar)');
    LEGEND_RATES(:,69) = strpad('d/dt conc_NO_cyt in component NO_Kinetics (Molar)');
    LEGEND_RATES(:,70) = strpad('d/dt conc_O2_cyt in component O2_Kinetics (Molar)');
    LEGEND_RATES(:,68) = strpad('d/dt conc_H2O2_cyt in component H2O2_Kinetics (Molar)');
    LEGEND_RATES(:,71) = strpad('d/dt conc_Eb in component sGC_Activation (dimensionless)');
    LEGEND_RATES(:,72) = strpad('d/dt conc_E5c in component sGC_Activation (dimensionless)');
    LEGEND_RATES(:,67) = strpad('d/dt conc_cGMP in component cGMP_Kinetics (Molar)');
    LEGEND_RATES(:,73) = strpad('d/dt conc_CaM___ in component MLCK_Activation_CaM_md (Molar)');
    LEGEND_RATES(:,19) = strpad('d/dt conc_CaMNCM in component MLCK_Activation_CaM_md (Molar)');
    LEGEND_STATES  = LEGEND_STATES';
    LEGEND_ALGEBRAIC = LEGEND_ALGEBRAIC';
    LEGEND_RATES = LEGEND_RATES';
    LEGEND_CONSTANTS = LEGEND_CONSTANTS';
end

function [STATES, CONSTANTS] = initConsts()
    VOI = 0; CONSTANTS = []; STATES = []; ALGEBRAIC = [];
    CONSTANTS(:,1) = 96487;
    CONSTANTS(:,2) = 5E-16;
    CONSTANTS(:,3) = 3E-18;
    CONSTANTS(:,4) = 7E-17;
    CONSTANTS(:,5) = 3.5E-16;
    CONSTANTS(:,6) = 1.7;
    CONSTANTS(:,7) = 2;
    CONSTANTS(:,8) = 12;
    CONSTANTS(:,9) = 3.7224E-4;
    CONSTANTS(:,10) = 6.3601E-4;
    CONSTANTS(:,11) = 5.71E-5;
    CONSTANTS(:,12) = 0.1;
    CONSTANTS(:,13) = 20;
    CONSTANTS(:,14) = 1.85;
    CONSTANTS(:,15) = 2.4E-4;
    CONSTANTS(:,16) = 0.5;
    CONSTANTS(:,17) = 1.1E-3;
    CONSTANTS(:,18) = 1;
    CONSTANTS(:,19) = 1.3E-4;
    CONSTANTS(:,20) = 2E2;
    CONSTANTS(:,21) = 1.445E-4;
    CONSTANTS(:,22) = 2E4;
    CONSTANTS(:,23) = 8.234E-5;
    CONSTANTS(:,24) = 7.5E-11;
    CONSTANTS(:,25) = 0.1;
    CONSTANTS(:,26) = 8.78E-2;
    CONSTANTS(:,27) = 6E-17;
    CONSTANTS(:,28) = 3.5E-18;
    CONSTANTS(:,29) = 6E-17;
    CONSTANTS(:,30) = 46;
    CONSTANTS(:,31) = 4.0;
    CONSTANTS(:,32) = 78;
    CONSTANTS(:,33) = 78;
    CONSTANTS(:,34) = 3.0E3;
    CONSTANTS(:,35) = 50;
    CONSTANTS(:,36) = 10;
    CONSTANTS(:,37) = 30;
    CONSTANTS(:,38) = 300;
    CONSTANTS(:,39) = 1.2E-8;
    CONSTANTS(:,40) = 1E-11;
    CONSTANTS(:,41) = 1E-11;
    CONSTANTS(:,42) = 1.2E-11;
    CONSTANTS(:,43) = 1E-11;
    CONSTANTS(:,44) = 3E-10;
    CONSTANTS(:,45) = 5.4;
    CONSTANTS(:,46) = 9.8325E-9;
    CONSTANTS(:,47) = 0.463;
    CONSTANTS(:,48) = 1.036E-9;
    CONSTANTS(:,49) = 5.00E-10;
    CONSTANTS(:,50) = 5E-4;
    CONSTANTS(:,51) = 1.15E-2;
    CONSTANTS(:,52) = 1.5;
    CONSTANTS(:,53) = 12;
    CONSTANTS(:,54) = 22;
    CONSTANTS(:,55) = 1.7097E-10;
    CONSTANTS(:,56) = 3.59E-3;
    CONSTANTS(:,57) = 1.3;
    CONSTANTS(:,58) = 12.29;
    CONSTANTS(:,59) = 87.5;
    CONSTANTS(:,60) = 0.27;
    CONSTANTS(:,61) = 0.35;
    CONSTANTS(:,62) = 1.25E-4;
    CONSTANTS(:,63) = 2.860E-9;
    CONSTANTS(:,64) = 0.0001;
    CONSTANTS(:,65) = 0.001;
    CONSTANTS(:,66) = 1.7E-4;
    CONSTANTS(:,67) = 1.4151E-9;
    CONSTANTS(:,68) = 4E-1;
    CONSTANTS(:,69) = 100E3;
    CONSTANTS(:,70) = 300;
    CONSTANTS(:,71) = 1.4E-1;
    CONSTANTS(:,72) = 1E5;
    CONSTANTS(:,73) = 8E4;
    CONSTANTS(:,74) = 0.8580;
    CONSTANTS(:,75) = 0.1420;
    CONSTANTS(:,76) = 8.3144621;
    CONSTANTS(:,77) = 298;
    CONSTANTS(:,78) = 1.21E-9;
    CONSTANTS(:,79) = 5E-7;
    CONSTANTS(:,80) = 0.001;
    CONSTANTS(:,81) = 3E-10;
    CONSTANTS(:,82) = 3E-10;
    STATES(:,1) = 140.55;
    STATES(:,2) = 131.37;
    STATES(:,3) = 9.1964;
    STATES(:,4) = 20.742;
    STATES(:,5) = 23.257;
    STATES(:,6) = 16.737;
    STATES(:,7) = 1.0533E-4;
    STATES(:,8) = 1.1023E-5;
    STATES(:,9) = 1.3568E-2;
    STATES(:,10) = 4.1154E-5;
    STATES(:,11) = 1.9451E-3;
    STATES(:,12) = 5.3950E-5;
    STATES(:,13) = 7.9487E-3;
    STATES(:,14) = 4.4109E-2;
    STATES(:,15) = 3.7182E-4;
    STATES(:,16) = 1.2766E-4;
    STATES(:,17) = 4.7661E-4;
    STATES(:,18) = 1.4863E-3;
    STATES(:,19) = 5.1368E-4;
    STATES(:,20) = 5.6334E-1;
    STATES(:,21) = 4.9091E-1;
    STATES(:,22) = 5.7849E-2;
    STATES(:,23) = 5.3239E-2;
    STATES(:,24) = -7.8960E-2;
    STATES(:,25) = -7.6961E-2;
    CONSTANTS(:,83) = 5.4;
    CONSTANTS(:,84) = 140;
    CONSTANTS(:,85) = 110;
    CONSTANTS(:,86) = 2.0;
    CONSTANTS(:,87) = 1;
    CONSTANTS(:,88) = 1;
    CONSTANTS(:,89) = -1;
    CONSTANTS(:,90) = 2;
    STATES(:,26) = 0.11187;
    STATES(:,27) = 0.89971;
    STATES(:,28) = 2.4140E-4;
    STATES(:,29) = 2.7808E-4;
    STATES(:,30) = 0.11631;
    STATES(:,31) = 0.028019;
    STATES(:,32) = 0.084779;
    STATES(:,33) = 0.072097;
    STATES(:,34) = 0.14878;
    STATES(:,35) = 0.12652;
    CONSTANTS(:,91) = 9.434E-4;
    CONSTANTS(:,92) = 1;
    CONSTANTS(:,93) = 2;
    CONSTANTS(:,94) = 0.125;
    STATES(:,36) = 0.0062437;
    STATES(:,37) = 0.000039459;
    STATES(:,38) = 0.00000024937;
    STATES(:,39) = 0.0056570;
    STATES(:,40) = 0.00010725;
    STATES(:,41) = 0.0000067783;
    STATES(:,42) = 2.1346E-2;
    STATES(:,43) = 4.7620E-4;
    STATES(:,44) = 1.0624E-5;
    STATES(:,45) = 1.9731E-2;
    STATES(:,46) = 1.3205E-3;
    STATES(:,47) = 2.94594E-4;
    CONSTANTS(:,95) = -3.8;
    CONSTANTS(:,96) = -10.0;
    CONSTANTS(:,97) = 224;
    STATES(:,48) = 0.0048876;
    STATES(:,49) = 0.0056051;
    STATES(:,50) = 0.0048876;
    STATES(:,51) = 0.0056051;
    STATES(:,52) = 0.94813;
    STATES(:,53) = 0.67538;
    STATES(:,54) = 0.94813;
    STATES(:,55) = 0.67538;
    STATES(:,56) = 1.0429E-5;
    STATES(:,57) = 1.5817E-5;
    STATES(:,58) = 0.87853;
    STATES(:,59) = 0.84170;
    STATES(:,60) = 3.7841E-5;
    STATES(:,61) = 4.9277E-5;
    STATES(:,62) = 1;
    STATES(:,63) = 1;
    CONSTANTS(:,98) = 1;
    CONSTANTS(:,99) = 7.7598E-13;
    CONSTANTS(:,100) = 1;
    CONSTANTS(:,101) = 7.7598E-13;
    CONSTANTS(:,102) = 2;
    CONSTANTS(:,103) = 7.7598E-13;
    CONSTANTS(:,104) = -1;
    CONSTANTS(:,105) = 7.7598E-13;
    CONSTANTS(:,106) = 1.2100E-11;
    CONSTANTS(:,107) = 1.9E7;
    CONSTANTS(:,108) = 0.01;
    CONSTANTS(:,109) = 1.6E6;
    CONSTANTS(:,110) = 3.4E4;
    CONSTANTS(:,111) = 1E-3;
    CONSTANTS(:,112) = 20E-3;
    CONSTANTS(:,113) = 2.8E-8;
    CONSTANTS(:,114) = 4.5E-3;
    CONSTANTS(:,115) = 6E-8;
    CONSTANTS(:,116) = 3E-4;
    CONSTANTS(:,117) = 6E-4;
    CONSTANTS(:,118) = 6E-4;
    CONSTANTS(:,119) = 2E6;
    CONSTANTS(:,120) = 100;
    CONSTANTS(:,121) = 0.1;
    CONSTANTS(:,122) = 3E3;
    CONSTANTS(:,123) = 98;
    CONSTANTS(:,124) = 1.09E-3;
    CONSTANTS(:,125) = 0.032;
    CONSTANTS(:,126) = 2.0E-3;
    CONSTANTS(:,127) = 5.5E-3;
    CONSTANTS(:,128) = 5.0E-4;
    CONSTANTS(:,129) = 1.0E-3;
    CONSTANTS(:,130) = 5.5E-4;
    CONSTANTS(:,131) = 66.9;
    CONSTANTS(:,132) = 2E-1;
    CONSTANTS(:,133) = 100;
    CONSTANTS(:,134) = -53.7;
    CONSTANTS(:,135) = 283.7;
    CONSTANTS(:,136) = 30.8E-3;
    CONSTANTS(:,137) = 2.0E-2;
    CONSTANTS(:,138) = 6.6E-2;
    CONSTANTS(:,139) = 1.6E-3;
    CONSTANTS(:,140) = 0.56;
    CONSTANTS(:,141) = 1E-3;
    CONSTANTS(:,142) = 4.0;
    CONSTANTS(:,143) = 2.8E6;
    CONSTANTS(:,144) = 6;
    CONSTANTS(:,145) = 100E6;
    CONSTANTS(:,146) = 800;
    CONSTANTS(:,147) = 2.8E6;
    CONSTANTS(:,148) = 6.0;
    CONSTANTS(:,149) = 100E6;
    CONSTANTS(:,150) = 800;
    CONSTANTS(:,151) = 1000E3;
    CONSTANTS(:,152) = 20;
    CONSTANTS(:,153) = 12.5E6;
    CONSTANTS(:,154) = 5.0;
    CONSTANTS(:,155) = 1000E3;
    CONSTANTS(:,156) = 1;
    CONSTANTS(:,157) = 4E3;
    CONSTANTS(:,158) = 0.4;
    CONSTANTS(:,159) = 1.8;
    CONSTANTS(:,160) = 0.1;
    CONSTANTS(:,161) = 0.045;
    CONSTANTS(:,162) = 1.00E-2;
    CONSTANTS(:,163) = 2.0E-3;
    CONSTANTS(:,164) = 30E-3;
    CONSTANTS(:,165) = 3.0995E-4;
    CONSTANTS(:,166) = 2.853E-12;
    CONSTANTS(:,167) = 100E-12;
    CONSTANTS(:,168) = 2.7225E-11;
    CONSTANTS(:,169) = 2.7225E-11;
    STATES(:,64) = 8.0265E-3;
    STATES(:,65) = 6.3288E-4;
    STATES(:,66) = 1.8089E-2;
    CONSTANTS(:,170) = 5.5E-3;
    STATES(:,67) = 9.0341E-3;
    STATES(:,68) = 5.0006E-5;
    STATES(:,69) = 9.2798E-5;
    STATES(:,70) = 8.3402E-7;
    STATES(:,71) = 0.27479;
    STATES(:,72) = 0.21715;
    CONSTANTS(:,171) = 1E-3;
    CONSTANTS(:,172) = 0.5263;
    STATES(:,73) = 7.3876E-3;
    CONSTANTS(:,173) = 3.9796;
    CONSTANTS(:,174) = 97.35;
    CONSTANTS(:,175) = 522.25;
    CONSTANTS(:,176) = 329.81;
    CONSTANTS(:,177) = ( CONSTANTS(:,81).*1.96000)./0.790000;
    CONSTANTS(:,178) = ( CONSTANTS(:,81).*1.33000)./0.790000;
    CONSTANTS(:,179) = ( CONSTANTS(:,76).*CONSTANTS(:,77))./CONSTANTS(:,1);
    CONSTANTS(:,180) =  (1.00000./7.00000).*(exp(CONSTANTS(:,84)./67.3000) - 1.00000);
    CONSTANTS(:,181) = ( CONSTANTS(:,80).*CONSTANTS(:,81))./CONSTANTS(:,79);
    CONSTANTS(:,182) = ( CONSTANTS(:,80).*CONSTANTS(:,82))./CONSTANTS(:,79);
    CONSTANTS(:,183) = ( CONSTANTS(:,80).*CONSTANTS(:,177))./CONSTANTS(:,79);
    CONSTANTS(:,184) = ( CONSTANTS(:,80).*CONSTANTS(:,178))./CONSTANTS(:,79);
    CONSTANTS(:,185) =  CONSTANTS(:,103).*CONSTANTS(:,181);
    CONSTANTS(:,186) =  CONSTANTS(:,105).*CONSTANTS(:,182);
    CONSTANTS(:,187) =  CONSTANTS(:,99).*CONSTANTS(:,183);
    CONSTANTS(:,188) =  CONSTANTS(:,101).*CONSTANTS(:,184);
    if (isempty(STATES)), warning('Initial values for states not set');, end
end

function [RATES, ALGEBRAIC] = computeRates(VOI, STATES, CONSTANTS)
    global algebraicVariableCount;
    statesSize = size(STATES);
    statesColumnCount = statesSize(2);
    if ( statesColumnCount == 1)
        STATES = STATES';
        ALGEBRAIC = zeros(1, algebraicVariableCount);
        utilOnes = 1;
    else
        statesRowCount = statesSize(1);
        ALGEBRAIC = zeros(statesRowCount, algebraicVariableCount);
        RATES = zeros(statesRowCount, statesColumnCount);
        utilOnes = ones(statesRowCount, 1);
    end
    RATES(:,9) =  CONSTANTS(:,69).*STATES(:,7).*(CONSTANTS(:,68) - STATES(:,9)) -  CONSTANTS(:,70).*STATES(:,9);
    RATES(:,28) =  CONSTANTS(:,14).*CONSTANTS(:,15).*((STATES(:,7)+ (1.00000 - CONSTANTS(:,16)).*CONSTANTS(:,17))./(STATES(:,7)+CONSTANTS(:,17))) -  CONSTANTS(:,18).*STATES(:,28);
    RATES(:,29) =  CONSTANTS(:,14).*CONSTANTS(:,15).*((STATES(:,15)+ (1.00000 - CONSTANTS(:,16)).*CONSTANTS(:,17))./(STATES(:,15)+CONSTANTS(:,17))) -  CONSTANTS(:,18).*STATES(:,29);
    RATES(:,30) =   - CONSTANTS(:,20).*( STATES(:,7).*STATES(:,30) -  CONSTANTS(:,21).*STATES(:,32)) -  CONSTANTS(:,22).*( STATES(:,7).*STATES(:,30) -  CONSTANTS(:,23).*STATES(:,34));
    RATES(:,31) =   - CONSTANTS(:,20).*( STATES(:,15).*STATES(:,31) -  CONSTANTS(:,21).*STATES(:,33)) -  CONSTANTS(:,22).*( STATES(:,15).*STATES(:,31) -  CONSTANTS(:,23).*STATES(:,35));
    RATES(:,13) = ((  - CONSTANTS(:,143).*power(STATES(:,7), 2.00000).*STATES(:,13)+ CONSTANTS(:,144).*STATES(:,10)) -  CONSTANTS(:,145).*power(STATES(:,7), 2.00000).*STATES(:,13))+ CONSTANTS(:,146).*STATES(:,8);
    RATES(:,69) = ( (CONSTANTS(:,78)./( 2.00000.*CONSTANTS(:,2))).*(CONSTANTS(:,113) -  CONSTANTS(:,116).*STATES(:,69)) -  CONSTANTS(:,107).*STATES(:,69).*STATES(:,70)) -  CONSTANTS(:,108).*STATES(:,69);
    RATES(:,68) = ( (CONSTANTS(:,78)./( 2.00000.*CONSTANTS(:,2))).*(CONSTANTS(:,115) -  CONSTANTS(:,118).*STATES(:,68))+ CONSTANTS(:,109).*CONSTANTS(:,111).*STATES(:,70)) -  CONSTANTS(:,110).*STATES(:,68).*CONSTANTS(:,112);
    RATES(:,67) =  CONSTANTS(:,124).*STATES(:,72) - ( CONSTANTS(:,125).*power(STATES(:,67), 2.00000))./(CONSTANTS(:,126)+STATES(:,67));
    RATES(:,73) = ((  - CONSTANTS(:,143).*power(STATES(:,15), 2.00000).*STATES(:,73)+ CONSTANTS(:,144).*STATES(:,17)) -  CONSTANTS(:,145).*power(STATES(:,15), 2.00000).*STATES(:,73))+ CONSTANTS(:,146).*STATES(:,16);
    RATES(:,14) =  CONSTANTS(:,69).*STATES(:,15).*(CONSTANTS(:,68) - STATES(:,14)) -  CONSTANTS(:,70).*STATES(:,14);
    ALGEBRAIC(:,1) = (1.00000+power(CONSTANTS(:,9)./STATES(:,7), 4.00000)+power(STATES(:,7)./CONSTANTS(:,10), 3.00000))./(1.00000+1.00000./CONSTANTS(:,11)+power(CONSTANTS(:,9)./STATES(:,7), 4.00000)+power(STATES(:,7)./CONSTANTS(:,10), 3.00000));
    RATES(:,26) = ( CONSTANTS(:,12).*(ALGEBRAIC(:,1) - STATES(:,26)))./ALGEBRAIC(:,1);
    ALGEBRAIC(:,2) = (1.00000+power(CONSTANTS(:,9)./STATES(:,15), 4.00000)+power(STATES(:,15)./CONSTANTS(:,10), 3.00000))./(1.00000+1.00000./CONSTANTS(:,11)+power(CONSTANTS(:,9)./STATES(:,15), 4.00000)+power(STATES(:,15)./CONSTANTS(:,10), 3.00000));
    RATES(:,27) = ( CONSTANTS(:,12).*(ALGEBRAIC(:,2) - STATES(:,27)))./ALGEBRAIC(:,2);
    ALGEBRAIC(:,9) = 1.00000./(1.00000+exp( - (STATES(:,24)+0.0239000)./0.00480000));
    RATES(:,56) = (ALGEBRAIC(:,9) - STATES(:,56))./CONSTANTS(:,64);
    ALGEBRAIC(:,10) = 1.00000./(1.00000+exp( - (STATES(:,25)+0.0239000)./0.00480000));
    RATES(:,57) = (ALGEBRAIC(:,10) - STATES(:,57))./CONSTANTS(:,64);
    ALGEBRAIC(:,11) = 1.00000./(1.00000+exp((STATES(:,24)+0.0661000)./0.00650000));
    RATES(:,58) = (ALGEBRAIC(:,11) - STATES(:,58))./CONSTANTS(:,65);
    ALGEBRAIC(:,12) = 1.00000./(1.00000+exp((STATES(:,25)+0.0661000)./0.00650000));
    RATES(:,59) = (ALGEBRAIC(:,12) - STATES(:,59))./CONSTANTS(:,65);
    ALGEBRAIC(:,18) = (  - CONSTANTS(:,1).*STATES(:,24))./( CONSTANTS(:,76).*CONSTANTS(:,77));
    RATES(:,70) = ((CONSTANTS(:,114) -  (CONSTANTS(:,78)./( 2.00000.*CONSTANTS(:,2))).*CONSTANTS(:,117).*STATES(:,70).*(ALGEBRAIC(:,18)./(1.00000 - exp( - ALGEBRAIC(:,18))))) -  CONSTANTS(:,107).*STATES(:,69).*STATES(:,70)) -  CONSTANTS(:,109).*CONSTANTS(:,111).*STATES(:,70);
    ALGEBRAIC(:,20) =  CONSTANTS(:,97).*exp(CONSTANTS(:,95)+ CONSTANTS(:,96).*STATES(:,24));
    [CONSTANTS, STATES, ALGEBRAIC] = rootfind_0(VOI, CONSTANTS, STATES, ALGEBRAIC);
    RATES(:,36) = ((  - CONSTANTS(:,34).*STATES(:,7).*(STATES(:,36) - ALGEBRAIC(:,3))+ CONSTANTS(:,35).*(STATES(:,37) - STATES(:,36))) -  CONSTANTS(:,36).*STATES(:,36))+ ALGEBRAIC(:,20).*STATES(:,39);
    RATES(:,37) = ((  - CONSTANTS(:,34).*STATES(:,7).*(STATES(:,37) - STATES(:,36))+ CONSTANTS(:,35).*(STATES(:,38) - STATES(:,37))) -  CONSTANTS(:,37).*STATES(:,37))+ ALGEBRAIC(:,20).*STATES(:,40);
    RATES(:,38) = (( CONSTANTS(:,34).*STATES(:,7).*STATES(:,37) -  CONSTANTS(:,35).*STATES(:,38)) -  CONSTANTS(:,38).*STATES(:,38))+ ALGEBRAIC(:,20).*STATES(:,41);
    RATES(:,39) =  CONSTANTS(:,36).*STATES(:,36) -  ALGEBRAIC(:,20).*STATES(:,39);
    RATES(:,40) =  CONSTANTS(:,37).*STATES(:,37) -  ALGEBRAIC(:,20).*STATES(:,40);
    RATES(:,41) =  CONSTANTS(:,38).*STATES(:,38) -  ALGEBRAIC(:,20).*STATES(:,41);
    ALGEBRAIC(:,21) =  CONSTANTS(:,97).*exp(CONSTANTS(:,95)+ CONSTANTS(:,96).*STATES(:,25));
    [CONSTANTS, STATES, ALGEBRAIC] = rootfind_1(VOI, CONSTANTS, STATES, ALGEBRAIC);
    RATES(:,42) = ((  - CONSTANTS(:,34).*STATES(:,15).*(STATES(:,42) - ALGEBRAIC(:,4))+ CONSTANTS(:,35).*(STATES(:,43) - STATES(:,42))) -  CONSTANTS(:,36).*STATES(:,42))+ ALGEBRAIC(:,21).*STATES(:,45);
    RATES(:,43) = ((  - CONSTANTS(:,34).*STATES(:,15).*(STATES(:,43) - STATES(:,42))+ CONSTANTS(:,35).*(STATES(:,44) - STATES(:,43))) -  CONSTANTS(:,37).*STATES(:,43))+ ALGEBRAIC(:,21).*STATES(:,46);
    RATES(:,44) = (( CONSTANTS(:,34).*STATES(:,15).*STATES(:,43) -  CONSTANTS(:,35).*STATES(:,44)) -  CONSTANTS(:,38).*STATES(:,44))+ ALGEBRAIC(:,21).*STATES(:,47);
    RATES(:,45) =  CONSTANTS(:,36).*STATES(:,42) -  ALGEBRAIC(:,21).*STATES(:,45);
    RATES(:,46) =  CONSTANTS(:,37).*STATES(:,43) -  ALGEBRAIC(:,21).*STATES(:,46);
    RATES(:,47) =  CONSTANTS(:,38).*STATES(:,44) -  ALGEBRAIC(:,21).*STATES(:,47);
    ALGEBRAIC(:,5) = 1.00000./(1.00000+exp( - (STATES(:,24)+0.00177000)./0.0145200));
    ALGEBRAIC(:,22) =  0.210987.*exp( - power((STATES(:,24)+0.214340)./0.195350, 2.00000)) - 0.0205900;
    RATES(:,48) = (ALGEBRAIC(:,5) - STATES(:,48))./ALGEBRAIC(:,22);
    ALGEBRAIC(:,6) = 1.00000./(1.00000+exp( - (STATES(:,25)+0.00177000)./0.0145200));
    ALGEBRAIC(:,23) =  0.210987.*exp( - power((STATES(:,25)+0.214340)./0.195350, 2.00000)) - 0.0205900;
    RATES(:,49) = (ALGEBRAIC(:,6) - STATES(:,49))./ALGEBRAIC(:,23);
    ALGEBRAIC(:,24) =  0.821390.*exp( - power((STATES(:,24)+0.0315900)./0.0274600, 2.00000))+0.000190000;
    RATES(:,50) = (ALGEBRAIC(:,5) - STATES(:,50))./ALGEBRAIC(:,24);
    ALGEBRAIC(:,25) =  0.821390.*exp( - power((STATES(:,25)+0.0315900)./0.0274600, 2.00000))+0.000190000;
    RATES(:,51) = (ALGEBRAIC(:,6) - STATES(:,51))./ALGEBRAIC(:,25);
    ALGEBRAIC(:,7) =  ((( CONSTANTS(:,134).*arbitrary_log(STATES(:,7)./1.00000, 10) - CONSTANTS(:,135)) -  CONSTANTS(:,133).*(STATES(:,69)./(CONSTANTS(:,132)+STATES(:,69)))) -  CONSTANTS(:,131).*(power(STATES(:,67), 2.00000)./(power(CONSTANTS(:,130), 2.00000)+power(STATES(:,67), 2.00000)))).*CONSTANTS(:,171);
    ALGEBRAIC(:,26) = 1.00000./(1.00000+exp( - (STATES(:,24) - ALGEBRAIC(:,7))./0.0308000));
    RATES(:,52) = (ALGEBRAIC(:,26) - STATES(:,52))./CONSTANTS(:,50);
    ALGEBRAIC(:,8) =  ((( CONSTANTS(:,134).*arbitrary_log(STATES(:,15)./1.00000, 10) - CONSTANTS(:,135)) -  CONSTANTS(:,133).*(STATES(:,69)./(CONSTANTS(:,132)+STATES(:,69)))) -  CONSTANTS(:,131).*(power(STATES(:,67), 2.00000)./(power(CONSTANTS(:,130), 2.00000)+power(STATES(:,67), 2.00000)))).*CONSTANTS(:,171);
    ALGEBRAIC(:,27) = 1.00000./(1.00000+exp( - (STATES(:,25) - ALGEBRAIC(:,8))./0.0308000));
    RATES(:,53) = (ALGEBRAIC(:,27) - STATES(:,53))./CONSTANTS(:,50);
    RATES(:,54) = (ALGEBRAIC(:,26) - STATES(:,54))./CONSTANTS(:,51);
    RATES(:,55) = (ALGEBRAIC(:,27) - STATES(:,55))./CONSTANTS(:,51);
    ALGEBRAIC(:,13) = 1.00000./(1.00000+exp( - (STATES(:,24)+0.00188000)./0.00757000));
    ALGEBRAIC(:,28) =  0.00289000.*exp( - power((STATES(:,24)+0.00863000)./0.0123900, 2.00000))+0.00243000;
    RATES(:,60) = (ALGEBRAIC(:,13) - STATES(:,60))./ALGEBRAIC(:,28);
    ALGEBRAIC(:,14) = 1.00000./(1.00000+exp( - (STATES(:,25)+0.00188000)./0.00757000));
    ALGEBRAIC(:,29) =  0.0289000.*exp( - power((STATES(:,25)+0.00863000)./0.0123900, 2.00000))+0.00243000;
    RATES(:,61) = (ALGEBRAIC(:,14) - STATES(:,61))./ALGEBRAIC(:,29);
    ALGEBRAIC(:,15) = 1.00000./(1.00000+exp((STATES(:,24)+0.0293200)./0.00154000));
    ALGEBRAIC(:,30) =  0.295590.*exp( - power((STATES(:,24) - 0.00472000)./0.112550, 2.00000))+0.0231900;
    RATES(:,62) = (ALGEBRAIC(:,15) - STATES(:,62))./ALGEBRAIC(:,30);
    ALGEBRAIC(:,16) = 1.00000./(1.00000+exp((STATES(:,25)+0.0293200)./0.00154000));
    ALGEBRAIC(:,31) =  0.295590.*exp( - power((STATES(:,25) - 0.00472000)./0.112550, 2.00000))+0.0231900;
    RATES(:,63) = (ALGEBRAIC(:,16) - STATES(:,63))./ALGEBRAIC(:,31);
    ALGEBRAIC(:,19) = 1.00000 - (STATES(:,71)+STATES(:,72));
    ALGEBRAIC(:,33) =  CONSTANTS(:,123).*STATES(:,67);
    RATES(:,71) =   - CONSTANTS(:,119).*STATES(:,71).*STATES(:,69)+ CONSTANTS(:,120).*ALGEBRAIC(:,19)+ ALGEBRAIC(:,33).*STATES(:,72);
    RATES(:,72) = ( CONSTANTS(:,122).*ALGEBRAIC(:,19).*STATES(:,69)+ CONSTANTS(:,121).*ALGEBRAIC(:,19)) -  ALGEBRAIC(:,33).*STATES(:,72);
    ALGEBRAIC(:,32) =  CONSTANTS(:,157).*STATES(:,12);
    ALGEBRAIC(:,36) = 1.00000+ CONSTANTS(:,142).*(power(STATES(:,68), 2.00000)./(power(CONSTANTS(:,141), 2.00000)+power(STATES(:,68), 2.00000)));
    ALGEBRAIC(:,39) = ( CONSTANTS(:,158).*(1.00000+ 3.65000.*(power(STATES(:,67), 2.00000)./(power(CONSTANTS(:,170), 2.00000)+power(STATES(:,67), 2.00000)))))./ALGEBRAIC(:,36);
    RATES(:,64) =   - ALGEBRAIC(:,32).*STATES(:,64)+ ALGEBRAIC(:,39).*STATES(:,65)+ CONSTANTS(:,161).*STATES(:,66);
    ALGEBRAIC(:,17) = CONSTANTS(:,164) - (STATES(:,64)+STATES(:,65)+STATES(:,66));
    RATES(:,65) = (( ALGEBRAIC(:,32).*STATES(:,64) -  ALGEBRAIC(:,39).*STATES(:,65)) -  CONSTANTS(:,159).*STATES(:,65))+ CONSTANTS(:,160).*ALGEBRAIC(:,17);
    ALGEBRAIC(:,41) = ALGEBRAIC(:,39);
    RATES(:,66) = (  - CONSTANTS(:,161).*STATES(:,66)+ ALGEBRAIC(:,41).*ALGEBRAIC(:,17)) -  ALGEBRAIC(:,32).*STATES(:,66);
    ALGEBRAIC(:,63) =  (STATES(:,28)./CONSTANTS(:,19)).*STATES(:,34);
    ALGEBRAIC(:,61) =  (STATES(:,28)./CONSTANTS(:,91)).*STATES(:,32);
    ALGEBRAIC(:,59) =  (STATES(:,28)./CONSTANTS(:,19)).*STATES(:,30);
    [CONSTANTS, STATES, ALGEBRAIC] = rootfind_2(VOI, CONSTANTS, STATES, ALGEBRAIC);
    RATES(:,32) =  CONSTANTS(:,20).*( STATES(:,7).*STATES(:,30) -  CONSTANTS(:,21).*STATES(:,32)) -  CONSTANTS(:,22).*( STATES(:,7).*STATES(:,32) -  CONSTANTS(:,23).*ALGEBRAIC(:,65));
    RATES(:,34) =  CONSTANTS(:,22).*( STATES(:,7).*STATES(:,30) -  CONSTANTS(:,23).*STATES(:,34)) -  CONSTANTS(:,20).*( STATES(:,7).*STATES(:,34) -  CONSTANTS(:,21).*ALGEBRAIC(:,65));
    ALGEBRAIC(:,67) =  (STATES(:,29)./CONSTANTS(:,19)).*STATES(:,35);
    ALGEBRAIC(:,62) =  (STATES(:,29)./CONSTANTS(:,91)).*STATES(:,33);
    ALGEBRAIC(:,60) =  (STATES(:,29)./CONSTANTS(:,19)).*STATES(:,31);
    [CONSTANTS, STATES, ALGEBRAIC] = rootfind_3(VOI, CONSTANTS, STATES, ALGEBRAIC);
    RATES(:,33) =  CONSTANTS(:,20).*( STATES(:,15).*STATES(:,31) -  CONSTANTS(:,21).*STATES(:,33)) -  CONSTANTS(:,22).*( STATES(:,15).*STATES(:,33) -  CONSTANTS(:,23).*ALGEBRAIC(:,69));
    RATES(:,35) =  CONSTANTS(:,22).*( STATES(:,15).*STATES(:,31) -  CONSTANTS(:,23).*STATES(:,35)) -  CONSTANTS(:,20).*( STATES(:,15).*STATES(:,35) -  CONSTANTS(:,21).*ALGEBRAIC(:,69));
    RATES(:,22) =  CONSTANTS(:,72).*STATES(:,20).*(CONSTANTS(:,71) - STATES(:,22)) -  CONSTANTS(:,73).*STATES(:,22);
    RATES(:,23) =  CONSTANTS(:,72).*STATES(:,21).*(CONSTANTS(:,71) - STATES(:,23)) -  CONSTANTS(:,73).*STATES(:,23);
    ALGEBRAIC(:,35) =  (CONSTANTS(:,179)./CONSTANTS(:,87)).*log(CONSTANTS(:,83)./STATES(:,1));
    ALGEBRAIC(:,47) =  (CONSTANTS(:,179)./CONSTANTS(:,89)).*log(CONSTANTS(:,85)./STATES(:,5));
    ALGEBRAIC(:,78) =  CONSTANTS(:,28).*CONSTANTS(:,74).*((ALGEBRAIC(:,35) - ALGEBRAIC(:,47))./((ALGEBRAIC(:,35) - ALGEBRAIC(:,47))+CONSTANTS(:,26)));
    ALGEBRAIC(:,80) =  CONSTANTS(:,29).*CONSTANTS(:,74).*(( ( CONSTANTS(:,84).*CONSTANTS(:,83).*power(CONSTANTS(:,85), 2.00000) -  STATES(:,3).*STATES(:,1).*power(STATES(:,5), 2.00000)).*CONSTANTS(:,30).*CONSTANTS(:,31).*CONSTANTS(:,32).*CONSTANTS(:,33))./( (CONSTANTS(:,84)+CONSTANTS(:,30)).*(CONSTANTS(:,83)+CONSTANTS(:,31)).*(CONSTANTS(:,85)+CONSTANTS(:,32)).*(CONSTANTS(:,85)+CONSTANTS(:,33)).*(STATES(:,3)+CONSTANTS(:,30)).*(STATES(:,1)+CONSTANTS(:,31)).*(STATES(:,5)+CONSTANTS(:,32)).*(STATES(:,5)+CONSTANTS(:,33))));
    ALGEBRAIC(:,43) =  (CONSTANTS(:,179)./CONSTANTS(:,88)).*log(CONSTANTS(:,84)./STATES(:,3));
    ALGEBRAIC(:,76) =  CONSTANTS(:,27).*CONSTANTS(:,74).*(power(ALGEBRAIC(:,43) - ALGEBRAIC(:,47), 4.00000)./(power(ALGEBRAIC(:,43) - ALGEBRAIC(:,47), 4.00000)+power(CONSTANTS(:,26), 4.00000)));
    ALGEBRAIC(:,48) =  CONSTANTS(:,43).*CONSTANTS(:,74).*(STATES(:,24) - ALGEBRAIC(:,47));
    ALGEBRAIC(:,82) =  CONSTANTS(:,39).*CONSTANTS(:,74).*(STATES(:,39)+STATES(:,40)+STATES(:,41)).*(STATES(:,24) - ALGEBRAIC(:,47));
    ALGEBRAIC(:,121) =  (( CONSTANTS(:,104).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77))).*(STATES(:,25) - STATES(:,24));
    ALGEBRAIC(:,122) = piecewise({STATES(:,25) - STATES(:,24)==0.00000,  CONSTANTS(:,186).*(STATES(:,6) - STATES(:,5)) },  CONSTANTS(:,186).*ALGEBRAIC(:,121).*((STATES(:,6) -  STATES(:,5).*exp( - ALGEBRAIC(:,121)))./(1.00000 - exp( - ALGEBRAIC(:,121)))));
    RATES(:,5) = (ALGEBRAIC(:,48)+ALGEBRAIC(:,82)+( CONSTANTS(:,1).*ALGEBRAIC(:,76)+ CONSTANTS(:,1).*(ALGEBRAIC(:,78)+ 2.00000.*ALGEBRAIC(:,80)+ALGEBRAIC(:,122))))./( CONSTANTS(:,1).*CONSTANTS(:,2));
    ALGEBRAIC(:,40) =  (CONSTANTS(:,179)./CONSTANTS(:,87)).*log(CONSTANTS(:,83)./STATES(:,2));
    ALGEBRAIC(:,49) =  (CONSTANTS(:,179)./CONSTANTS(:,89)).*log(CONSTANTS(:,85)./STATES(:,6));
    ALGEBRAIC(:,79) =  CONSTANTS(:,28).*CONSTANTS(:,75).*((ALGEBRAIC(:,40) - ALGEBRAIC(:,49))./((ALGEBRAIC(:,40) - ALGEBRAIC(:,49))+CONSTANTS(:,26)));
    ALGEBRAIC(:,81) =  CONSTANTS(:,29).*CONSTANTS(:,75).*(( ( CONSTANTS(:,84).*CONSTANTS(:,83).*power(CONSTANTS(:,85), 2.00000) -  STATES(:,4).*STATES(:,2).*power(STATES(:,6), 2.00000)).*CONSTANTS(:,30).*CONSTANTS(:,31).*CONSTANTS(:,32).*CONSTANTS(:,33))./( (CONSTANTS(:,84)+CONSTANTS(:,30)).*(CONSTANTS(:,83)+CONSTANTS(:,31)).*(CONSTANTS(:,85)+CONSTANTS(:,32)).*(CONSTANTS(:,85)+CONSTANTS(:,33)).*(STATES(:,4)+CONSTANTS(:,30)).*(STATES(:,2)+CONSTANTS(:,31)).*(STATES(:,6)+CONSTANTS(:,32)).*(STATES(:,6)+CONSTANTS(:,33))));
    ALGEBRAIC(:,45) =  (CONSTANTS(:,179)./CONSTANTS(:,88)).*log(CONSTANTS(:,84)./STATES(:,4));
    ALGEBRAIC(:,77) =  CONSTANTS(:,27).*CONSTANTS(:,75).*(power(ALGEBRAIC(:,45) - ALGEBRAIC(:,49), 4.00000)./(power(ALGEBRAIC(:,45) - ALGEBRAIC(:,49), 4.00000)+power(CONSTANTS(:,26), 4.00000)));
    ALGEBRAIC(:,50) =  CONSTANTS(:,43).*CONSTANTS(:,75).*(STATES(:,25) - ALGEBRAIC(:,49));
    ALGEBRAIC(:,83) =  CONSTANTS(:,39).*CONSTANTS(:,75).*(STATES(:,45)+STATES(:,46)+STATES(:,47)).*(STATES(:,25) - ALGEBRAIC(:,49));
    RATES(:,6) = (ALGEBRAIC(:,50)+ALGEBRAIC(:,83)+( CONSTANTS(:,1).*ALGEBRAIC(:,77)+ CONSTANTS(:,1).*((ALGEBRAIC(:,79)+ 2.00000.*ALGEBRAIC(:,81)) - ALGEBRAIC(:,122))))./( CONSTANTS(:,1).*CONSTANTS(:,3));
    ALGEBRAIC(:,123) = CONSTANTS(:,162) - (STATES(:,13)+STATES(:,8)+STATES(:,10)+STATES(:,11)+STATES(:,12));
    RATES(:,8) = (( CONSTANTS(:,145).*power(STATES(:,7), 2.00000).*STATES(:,13) -  CONSTANTS(:,146).*STATES(:,8)) -  CONSTANTS(:,147).*power(STATES(:,7), 2.00000).*STATES(:,8))+ CONSTANTS(:,148).*ALGEBRAIC(:,123);
    ALGEBRAIC(:,124) = CONSTANTS(:,163) - (STATES(:,11)+STATES(:,12));
    RATES(:,10) = (((( CONSTANTS(:,143).*power(STATES(:,7), 2.00000).*STATES(:,13) -  CONSTANTS(:,144).*STATES(:,10)) -  CONSTANTS(:,149).*power(STATES(:,7), 2.00000).*STATES(:,10))+ CONSTANTS(:,150).*ALGEBRAIC(:,123)) -  CONSTANTS(:,151).*ALGEBRAIC(:,124).*STATES(:,10))+ CONSTANTS(:,152).*STATES(:,11);
    RATES(:,12) = (( CONSTANTS(:,153).*power(STATES(:,7), 2.00000).*STATES(:,11) -  CONSTANTS(:,154).*STATES(:,12))+ CONSTANTS(:,155).*ALGEBRAIC(:,124).*ALGEBRAIC(:,123)) -  CONSTANTS(:,156).*STATES(:,12);
    RATES(:,11) = (( CONSTANTS(:,151).*ALGEBRAIC(:,124).*STATES(:,10) -  CONSTANTS(:,152).*STATES(:,11)) -  CONSTANTS(:,153).*power(STATES(:,7), 2.00000).*STATES(:,11))+ CONSTANTS(:,154).*STATES(:,12);
    ALGEBRAIC(:,129) = power(STATES(:,70), 1.0 ./ 2)./(power(CONSTANTS(:,137), 1.0 ./ 2)+power(STATES(:,70), 1.0 ./ 2));
    ALGEBRAIC(:,130) = STATES(:,68)./(CONSTANTS(:,138)+STATES(:,68));
    ALGEBRAIC(:,131) =  CONSTANTS(:,167).*(1.00000 - ALGEBRAIC(:,129)).*(1.00000 - ALGEBRAIC(:,130));
    ALGEBRAIC(:,125) =  (CONSTANTS(:,165)./CONSTANTS(:,172)).*(1.00000 - STATES(:,67)./( 2.00000.*(CONSTANTS(:,128)+STATES(:,67))));
    ALGEBRAIC(:,132) =  ALGEBRAIC(:,131).*CONSTANTS(:,74).*((power(STATES(:,7)./ALGEBRAIC(:,125), CONSTANTS(:,7)) - power(STATES(:,20)./CONSTANTS(:,6), CONSTANTS(:,7)))./(1.00000+(power(STATES(:,7)./ALGEBRAIC(:,125), CONSTANTS(:,7))+power(STATES(:,20)./CONSTANTS(:,6), CONSTANTS(:,7)))));
    ALGEBRAIC(:,64) =  CONSTANTS(:,13).*power(ALGEBRAIC(:,63), 3.00000).*(STATES(:,20) - STATES(:,7)).*2.00000.*CONSTANTS(:,1).*CONSTANTS(:,4);
    ALGEBRAIC(:,68) =  CONSTANTS(:,13).*power(ALGEBRAIC(:,67), 3.00000).*(STATES(:,20) - STATES(:,15)).*2.00000.*CONSTANTS(:,1).*CONSTANTS(:,4);
    ALGEBRAIC(:,133) =  ALGEBRAIC(:,131).*CONSTANTS(:,75).*((power(STATES(:,15)./ALGEBRAIC(:,125), CONSTANTS(:,7)) - power(STATES(:,20)./CONSTANTS(:,6), CONSTANTS(:,7)))./(1.00000+(power(STATES(:,15)./ALGEBRAIC(:,125), CONSTANTS(:,7))+power(STATES(:,20)./CONSTANTS(:,6), CONSTANTS(:,7)))));
    RATES(:,20) = (((ALGEBRAIC(:,132) - ALGEBRAIC(:,64))+ALGEBRAIC(:,133)) - ALGEBRAIC(:,68))./( CONSTANTS(:,1).*CONSTANTS(:,4)) - RATES(:,22);
    ALGEBRAIC(:,51) =  (CONSTANTS(:,179)./CONSTANTS(:,90)).*log(CONSTANTS(:,86)./STATES(:,7));
    ALGEBRAIC(:,52) =  CONSTANTS(:,42).*CONSTANTS(:,74).*(STATES(:,24) - ALGEBRAIC(:,51));
    ALGEBRAIC(:,126) =  CONSTANTS(:,166).*(1.00000+( 3.60000.*STATES(:,67))./(CONSTANTS(:,129)+STATES(:,67)));
    ALGEBRAIC(:,127) =  ALGEBRAIC(:,126).*CONSTANTS(:,74).*(STATES(:,7)./(CONSTANTS(:,66)+STATES(:,7)));
    ALGEBRAIC(:,71) = (STATES(:,20)+STATES(:,21))./2.00000;
    ALGEBRAIC(:,72) =  CONSTANTS(:,24).*(CONSTANTS(:,25)./(CONSTANTS(:,25)+ALGEBRAIC(:,71))).*(STATES(:,24) - ALGEBRAIC(:,51));
    ALGEBRAIC(:,111) =  0.740000.*STATES(:,62)+0.260000;
    ALGEBRAIC(:,112) =  CONSTANTS(:,67).*CONSTANTS(:,74).*STATES(:,60).*ALGEBRAIC(:,111).*(STATES(:,24) - ALGEBRAIC(:,51));
    ALGEBRAIC(:,55) = ( STATES(:,26).*(1.00000+power(STATES(:,7)./CONSTANTS(:,10), 3.00000)))./(1.00000+power(CONSTANTS(:,9)./STATES(:,7), 4.00000)+power(STATES(:,7)./CONSTANTS(:,10), 3.00000));
    ALGEBRAIC(:,56) =  CONSTANTS(:,8).*ALGEBRAIC(:,55).*(STATES(:,21) - STATES(:,7)).*( 2.00000.*CONSTANTS(:,1).*CONSTANTS(:,4));
    ALGEBRAIC(:,119) =  (( CONSTANTS(:,102).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77))).*(STATES(:,25) - STATES(:,24));
    ALGEBRAIC(:,120) = piecewise({STATES(:,25) - STATES(:,24)==0.00000,  CONSTANTS(:,185).*(STATES(:,15) - STATES(:,7)) },  CONSTANTS(:,185).*ALGEBRAIC(:,119).*((STATES(:,15) -  STATES(:,7).*exp( - ALGEBRAIC(:,119)))./(1.00000 - exp( - ALGEBRAIC(:,119)))));
    ALGEBRAIC(:,134) =  ALGEBRAIC(:,131).*CONSTANTS(:,74).*((power(STATES(:,7)./ALGEBRAIC(:,125), CONSTANTS(:,7)) - power(STATES(:,21)./CONSTANTS(:,6), CONSTANTS(:,7)))./(1.00000+(power(STATES(:,7)./ALGEBRAIC(:,125), CONSTANTS(:,7))+power(STATES(:,21)./CONSTANTS(:,6), CONSTANTS(:,7)))));
    RATES(:,7) = (( - (((ALGEBRAIC(:,52)+ALGEBRAIC(:,127)+ALGEBRAIC(:,72)+ALGEBRAIC(:,112)+ALGEBRAIC(:,132)+ALGEBRAIC(:,134)) - ALGEBRAIC(:,56)) - ALGEBRAIC(:,64))+ 2.00000.*CONSTANTS(:,1).*ALGEBRAIC(:,120))./( 2.00000.*CONSTANTS(:,1).*CONSTANTS(:,5)) -  2.00000.*(RATES(:,8)+RATES(:,10)+RATES(:,11))) - RATES(:,9);
    ALGEBRAIC(:,135) =  ALGEBRAIC(:,131).*CONSTANTS(:,75).*((power(STATES(:,15)./ALGEBRAIC(:,125), CONSTANTS(:,7)) - power(STATES(:,21)./CONSTANTS(:,6), CONSTANTS(:,7)))./(1.00000+(power(STATES(:,15)./ALGEBRAIC(:,125), CONSTANTS(:,7))+power(STATES(:,21)./CONSTANTS(:,6), CONSTANTS(:,7)))));
    ALGEBRAIC(:,57) = ( STATES(:,27).*(1.00000+power(STATES(:,15)./CONSTANTS(:,10), 3.00000)))./(1.00000+power(CONSTANTS(:,9)./STATES(:,15), 4.00000)+power(STATES(:,15)./CONSTANTS(:,10), 3.00000));
    ALGEBRAIC(:,58) =  CONSTANTS(:,8).*ALGEBRAIC(:,57).*(STATES(:,21) - STATES(:,15)).*( 2.00000.*CONSTANTS(:,1).*CONSTANTS(:,4));
    RATES(:,21) = (((ALGEBRAIC(:,134) - ALGEBRAIC(:,56))+ALGEBRAIC(:,135)) - ALGEBRAIC(:,58))./( CONSTANTS(:,1).*CONSTANTS(:,4)) - RATES(:,23);
    ALGEBRAIC(:,137) = CONSTANTS(:,162) - (STATES(:,73)+STATES(:,16)+STATES(:,17)+STATES(:,18)+STATES(:,19));
    RATES(:,16) = (( CONSTANTS(:,145).*power(STATES(:,15), 2.00000).*STATES(:,73) -  CONSTANTS(:,146).*STATES(:,16)) -  CONSTANTS(:,147).*power(STATES(:,15), 2.00000).*STATES(:,16))+ CONSTANTS(:,148).*ALGEBRAIC(:,137);
    ALGEBRAIC(:,140) = CONSTANTS(:,163) - (STATES(:,18)+STATES(:,19));
    RATES(:,17) = (((( CONSTANTS(:,143).*power(STATES(:,15), 2.00000).*STATES(:,73) -  CONSTANTS(:,144).*STATES(:,17)) -  CONSTANTS(:,149).*power(STATES(:,15), 2.00000).*STATES(:,17))+ CONSTANTS(:,150).*ALGEBRAIC(:,137)) -  CONSTANTS(:,151).*ALGEBRAIC(:,140).*STATES(:,17))+ CONSTANTS(:,152).*STATES(:,18);
    RATES(:,19) = (( CONSTANTS(:,153).*power(STATES(:,15), 2.00000).*STATES(:,18) -  CONSTANTS(:,154).*STATES(:,19))+ CONSTANTS(:,155).*ALGEBRAIC(:,140).*ALGEBRAIC(:,137)) -  CONSTANTS(:,156).*STATES(:,19);
    ALGEBRAIC(:,38) =  CONSTANTS(:,40).*CONSTANTS(:,74).*(STATES(:,24) - ALGEBRAIC(:,35));
    ALGEBRAIC(:,84) = 3.24100./(1.00000+exp(((STATES(:,24) - ALGEBRAIC(:,35)) - 185.720)./1.59600));
    ALGEBRAIC(:,86) =  1.00000.*(( 0.0130000.*exp(((STATES(:,24) - ALGEBRAIC(:,35))+5.79220)./0.294000)+exp(((STATES(:,24) - ALGEBRAIC(:,35)) - 653.733)./0.244000))./(1.00000+exp( - ((STATES(:,24) - ALGEBRAIC(:,35)) - 7.77750)./0.0928000)));
    ALGEBRAIC(:,87) = ALGEBRAIC(:,84)./(ALGEBRAIC(:,84)+ALGEBRAIC(:,86));
    ALGEBRAIC(:,88) =  CONSTANTS(:,44).*CONSTANTS(:,74).*ALGEBRAIC(:,87).*power((CONSTANTS(:,83)./CONSTANTS(:,45)), 1.0 ./ 2).*(STATES(:,24) - ALGEBRAIC(:,35));
    ALGEBRAIC(:,96) =  CONSTANTS(:,48).*CONSTANTS(:,74).*power(CONSTANTS(:,83)./CONSTANTS(:,45), CONSTANTS(:,47)).*(STATES(:,24) - ALGEBRAIC(:,35));
    ALGEBRAIC(:,92) =  0.580000.*STATES(:,48)+ 0.420000.*STATES(:,50);
    ALGEBRAIC(:,93) =  CONSTANTS(:,46).*CONSTANTS(:,74).*power(ALGEBRAIC(:,92), 2.00000).*(STATES(:,24) - ALGEBRAIC(:,35));
    ALGEBRAIC(:,98) =  0.650000.*STATES(:,52)+ 0.350000.*STATES(:,54);
    ALGEBRAIC(:,99) =  CONSTANTS(:,49).*CONSTANTS(:,74).*ALGEBRAIC(:,98).*(STATES(:,24) - ALGEBRAIC(:,35));
    ALGEBRAIC(:,139) =  CONSTANTS(:,168).*(1.00000 - STATES(:,70)./(CONSTANTS(:,139)+STATES(:,70))).*(1.00000 - STATES(:,68)./(CONSTANTS(:,140)+STATES(:,68)));
    ALGEBRAIC(:,102) = 1.00000./(1.00000+ 0.124500.*exp((  - 0.100000.*STATES(:,24).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77)))+ 0.0365000.*CONSTANTS(:,180).*exp((  - STATES(:,24).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77))));
    ALGEBRAIC(:,141) =  ALGEBRAIC(:,139).*ALGEBRAIC(:,102).*(CONSTANTS(:,83)./(CONSTANTS(:,83)+CONSTANTS(:,52))).*(power(STATES(:,3), 1.50000)./(power(STATES(:,3), 1.50000)+power(CONSTANTS(:,53), 1.50000)));
    ALGEBRAIC(:,115) =  (( CONSTANTS(:,98).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77))).*(STATES(:,25) - STATES(:,24));
    ALGEBRAIC(:,116) = piecewise({STATES(:,25) - STATES(:,24)==0.00000,  CONSTANTS(:,187).*(STATES(:,2) - STATES(:,1)) },  CONSTANTS(:,187).*ALGEBRAIC(:,115).*((STATES(:,2) -  STATES(:,1).*exp( - ALGEBRAIC(:,115)))./(1.00000 - exp( - ALGEBRAIC(:,115)))));
    RATES(:,1) = ( - ((ALGEBRAIC(:,38)+ALGEBRAIC(:,88)+ALGEBRAIC(:,96)+ALGEBRAIC(:,93)+ALGEBRAIC(:,99)) -  2.00000.*ALGEBRAIC(:,141))+ CONSTANTS(:,1).*(ALGEBRAIC(:,78)+ALGEBRAIC(:,80)+ALGEBRAIC(:,116)))./( CONSTANTS(:,1).*CONSTANTS(:,2));
    ALGEBRAIC(:,42) =  CONSTANTS(:,40).*CONSTANTS(:,75).*(STATES(:,25) - ALGEBRAIC(:,40));
    ALGEBRAIC(:,85) = 3.24100./(1.00000+exp(((STATES(:,25) - ALGEBRAIC(:,40)) - 185.720)./1.59600));
    ALGEBRAIC(:,89) =  1.00000.*(( 0.0130000.*exp(((STATES(:,25) - ALGEBRAIC(:,40))+5.79220)./0.294000)+exp(((STATES(:,25) - ALGEBRAIC(:,40)) - 653.733)./0.244000))./(1.00000+exp( - (((STATES(:,25) - ALGEBRAIC(:,40)) - 7.77750)./0.0928000))));
    ALGEBRAIC(:,90) = ALGEBRAIC(:,85)./(ALGEBRAIC(:,85)+ALGEBRAIC(:,89));
    ALGEBRAIC(:,91) =  CONSTANTS(:,44).*CONSTANTS(:,75).*ALGEBRAIC(:,90).*power((CONSTANTS(:,83)./CONSTANTS(:,45)), 1.0 ./ 2).*(STATES(:,25) - ALGEBRAIC(:,40));
    ALGEBRAIC(:,97) =  CONSTANTS(:,48).*CONSTANTS(:,75).*power(CONSTANTS(:,83)./CONSTANTS(:,45), CONSTANTS(:,47)).*(STATES(:,25) - ALGEBRAIC(:,40));
    ALGEBRAIC(:,94) =  0.580000.*STATES(:,49)+ 0.420000.*STATES(:,51);
    ALGEBRAIC(:,95) =  CONSTANTS(:,46).*CONSTANTS(:,75).*power(ALGEBRAIC(:,94), 2.00000).*(STATES(:,25) - ALGEBRAIC(:,40));
    ALGEBRAIC(:,100) =  0.650000.*STATES(:,53)+ 0.350000.*STATES(:,55);
    ALGEBRAIC(:,101) =  CONSTANTS(:,49).*CONSTANTS(:,75).*ALGEBRAIC(:,100).*(STATES(:,25) - ALGEBRAIC(:,40));
    ALGEBRAIC(:,103) = 1.00000./(1.00000+ 0.124500.*exp((  - 0.100000.*STATES(:,25).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77)))+ 0.0365000.*CONSTANTS(:,180).*exp((  - STATES(:,25).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77))));
    ALGEBRAIC(:,142) =  ALGEBRAIC(:,139).*ALGEBRAIC(:,103).*(CONSTANTS(:,83)./(CONSTANTS(:,83)+CONSTANTS(:,52))).*(power(STATES(:,4), 1.50000)./(power(STATES(:,4), 1.50000)+power(CONSTANTS(:,54), 1.50000)));
    RATES(:,2) = ( - ((ALGEBRAIC(:,42)+ALGEBRAIC(:,91)+ALGEBRAIC(:,97)+ALGEBRAIC(:,95)+ALGEBRAIC(:,101)) -  2.00000.*ALGEBRAIC(:,142))+ CONSTANTS(:,1).*((ALGEBRAIC(:,79)+ALGEBRAIC(:,81)) - ALGEBRAIC(:,116)))./( CONSTANTS(:,1).*CONSTANTS(:,3));
    ALGEBRAIC(:,44) =  CONSTANTS(:,41).*CONSTANTS(:,74).*(STATES(:,24) - ALGEBRAIC(:,43));
    ALGEBRAIC(:,109) =  CONSTANTS(:,63).*CONSTANTS(:,74).*STATES(:,56).*STATES(:,58).*(STATES(:,24) - ALGEBRAIC(:,43));
    ALGEBRAIC(:,74) =  ALGEBRAIC(:,72).*( (power(CONSTANTS(:,92), 2.00000)./power(CONSTANTS(:,93), 2.00000)).*CONSTANTS(:,94)).*((STATES(:,3) -  CONSTANTS(:,84).*exp((  - CONSTANTS(:,92).*CONSTANTS(:,1).*STATES(:,24))./( CONSTANTS(:,76).*CONSTANTS(:,77))))./(STATES(:,7) -  CONSTANTS(:,86).*exp((  - CONSTANTS(:,93).*CONSTANTS(:,1).*STATES(:,24))./( CONSTANTS(:,76).*CONSTANTS(:,77))))).*((1.00000 - exp((  - CONSTANTS(:,93).*CONSTANTS(:,1).*STATES(:,24))./( CONSTANTS(:,76).*CONSTANTS(:,77))))./(1.00000 - exp((  - CONSTANTS(:,92).*CONSTANTS(:,1).*STATES(:,24))./( CONSTANTS(:,76).*CONSTANTS(:,77)))));
    ALGEBRAIC(:,117) =  (( CONSTANTS(:,100).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77))).*(STATES(:,25) - STATES(:,24));
    ALGEBRAIC(:,118) = piecewise({STATES(:,25) - STATES(:,24)==0.00000,  CONSTANTS(:,188).*(STATES(:,4) - STATES(:,3)) },  CONSTANTS(:,188).*ALGEBRAIC(:,117).*((STATES(:,4) -  STATES(:,3).*exp( - ALGEBRAIC(:,117)))./(1.00000 - exp( - ALGEBRAIC(:,117)))));
    RATES(:,3) = ( - (ALGEBRAIC(:,44)+ 3.00000.*ALGEBRAIC(:,141)+ALGEBRAIC(:,109)+ALGEBRAIC(:,74))+ CONSTANTS(:,1).*(ALGEBRAIC(:,76)+ALGEBRAIC(:,80)+ALGEBRAIC(:,118)))./( CONSTANTS(:,1).*CONSTANTS(:,2));
    ALGEBRAIC(:,46) =  CONSTANTS(:,41).*CONSTANTS(:,75).*(STATES(:,25) - ALGEBRAIC(:,45));
    ALGEBRAIC(:,110) =  CONSTANTS(:,63).*CONSTANTS(:,75).*STATES(:,57).*STATES(:,59).*(STATES(:,25) - ALGEBRAIC(:,45));
    ALGEBRAIC(:,53) =  (CONSTANTS(:,179)./CONSTANTS(:,90)).*log(CONSTANTS(:,86)./STATES(:,15));
    ALGEBRAIC(:,73) =  CONSTANTS(:,24).*(CONSTANTS(:,25)./(CONSTANTS(:,25)+ALGEBRAIC(:,71))).*(STATES(:,25) - ALGEBRAIC(:,53));
    ALGEBRAIC(:,75) =  ALGEBRAIC(:,73).*( (power(CONSTANTS(:,92), 2.00000)./power(CONSTANTS(:,93), 2.00000)).*CONSTANTS(:,94)).*((STATES(:,4) -  CONSTANTS(:,84).*exp((  - CONSTANTS(:,92).*CONSTANTS(:,1).*STATES(:,25))./( CONSTANTS(:,76).*CONSTANTS(:,77))))./(STATES(:,15) -  CONSTANTS(:,86).*exp((  - CONSTANTS(:,93).*CONSTANTS(:,1).*STATES(:,25))./( CONSTANTS(:,76).*CONSTANTS(:,77))))).*((1.00000 - exp((  - CONSTANTS(:,93).*CONSTANTS(:,1).*STATES(:,25))./( CONSTANTS(:,76).*CONSTANTS(:,77))))./(1.00000 - exp((  - CONSTANTS(:,92).*CONSTANTS(:,1).*STATES(:,25))./( CONSTANTS(:,76).*CONSTANTS(:,77)))));
    ALGEBRAIC(:,104) = 1.00000./(1.00000+power(CONSTANTS(:,62)./STATES(:,15), 2.00000));
    ALGEBRAIC(:,105) = exp(( CONSTANTS(:,61).*STATES(:,25).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77)));
    ALGEBRAIC(:,106) = exp(( (CONSTANTS(:,61) - 1.00000).*STATES(:,25).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77)));
    ALGEBRAIC(:,107) =  power(CONSTANTS(:,84), 3.00000).*STATES(:,15)+ power(STATES(:,4), 3.00000).*CONSTANTS(:,86)+ power(CONSTANTS(:,59), 3.00000).*STATES(:,15)+ CONSTANTS(:,57).*power(STATES(:,4), 3.00000)+ power(CONSTANTS(:,58), 3.00000).*CONSTANTS(:,86).*(1.00000+STATES(:,15)./CONSTANTS(:,56))+ CONSTANTS(:,56).*power(CONSTANTS(:,84), 3.00000).*(1.00000+power(STATES(:,4)./CONSTANTS(:,58), 3.00000));
    ALGEBRAIC(:,108) =  CONSTANTS(:,55).*ALGEBRAIC(:,104).*(( ALGEBRAIC(:,105).*power(STATES(:,4), 3.00000).*CONSTANTS(:,86) -  ALGEBRAIC(:,106).*power(CONSTANTS(:,84), 3.00000).*STATES(:,15))./( ALGEBRAIC(:,107).*(1.00000+ CONSTANTS(:,60).*ALGEBRAIC(:,106))));
    RATES(:,4) = ( - (ALGEBRAIC(:,46)+ 3.00000.*ALGEBRAIC(:,142)+ALGEBRAIC(:,110)+ALGEBRAIC(:,75)+ 3.00000.*ALGEBRAIC(:,108))+ CONSTANTS(:,1).*((ALGEBRAIC(:,77)+ALGEBRAIC(:,81)) - ALGEBRAIC(:,118)))./( CONSTANTS(:,1).*CONSTANTS(:,3));
    RATES(:,18) = (( CONSTANTS(:,151).*ALGEBRAIC(:,140).*STATES(:,17) -  CONSTANTS(:,152).*STATES(:,18)) -  CONSTANTS(:,153).*power(STATES(:,15), 2.00000).*STATES(:,18))+ CONSTANTS(:,154).*STATES(:,19);
    ALGEBRAIC(:,54) =  CONSTANTS(:,42).*CONSTANTS(:,75).*(STATES(:,25) - ALGEBRAIC(:,53));
    ALGEBRAIC(:,128) =  ALGEBRAIC(:,126).*CONSTANTS(:,75).*(STATES(:,15)./(CONSTANTS(:,66)+STATES(:,15)));
    ALGEBRAIC(:,113) =  0.740000.*STATES(:,63)+0.260000;
    ALGEBRAIC(:,114) =  CONSTANTS(:,67).*CONSTANTS(:,75).*STATES(:,61).*ALGEBRAIC(:,113).*(STATES(:,25) - ALGEBRAIC(:,53));
    RATES(:,15) = (( - (((((ALGEBRAIC(:,54)+ALGEBRAIC(:,128)+ALGEBRAIC(:,73)+ALGEBRAIC(:,114)) -  2.00000.*ALGEBRAIC(:,108))+ALGEBRAIC(:,135)+ALGEBRAIC(:,133)) - ALGEBRAIC(:,58)) - ALGEBRAIC(:,68)) -  2.00000.*CONSTANTS(:,1).*ALGEBRAIC(:,120))./( 2.00000.*CONSTANTS(:,1).*CONSTANTS(:,3)) -  2.00000.*(RATES(:,16)+RATES(:,17)+RATES(:,18))) - RATES(:,14);
    ALGEBRAIC(:,136) = ALGEBRAIC(:,132)+ALGEBRAIC(:,134);
    ALGEBRAIC(:,144) = ((((((((ALGEBRAIC(:,38)+ALGEBRAIC(:,88)+ALGEBRAIC(:,96)+ALGEBRAIC(:,93)+ALGEBRAIC(:,99)) -  CONSTANTS(:,1).*ALGEBRAIC(:,116))+ALGEBRAIC(:,141)+ALGEBRAIC(:,44)+ALGEBRAIC(:,109)+ALGEBRAIC(:,74)) -  CONSTANTS(:,1).*ALGEBRAIC(:,118))+ALGEBRAIC(:,48)+ALGEBRAIC(:,82)+ CONSTANTS(:,1).*ALGEBRAIC(:,122)+ALGEBRAIC(:,52)+ALGEBRAIC(:,127)+ALGEBRAIC(:,112)+ALGEBRAIC(:,72)) -  2.00000.*CONSTANTS(:,1).*ALGEBRAIC(:,120))+ALGEBRAIC(:,136)) - ALGEBRAIC(:,56)) - ALGEBRAIC(:,64);
    RATES(:,24) =  - ALGEBRAIC(:,144)./( CONSTANTS(:,74).*CONSTANTS(:,106));
    ALGEBRAIC(:,138) = ALGEBRAIC(:,133)+ALGEBRAIC(:,135);
    ALGEBRAIC(:,145) = ((((ALGEBRAIC(:,42)+ALGEBRAIC(:,91)+ALGEBRAIC(:,97)+ALGEBRAIC(:,95)+ALGEBRAIC(:,101)+ CONSTANTS(:,1).*ALGEBRAIC(:,116)+ALGEBRAIC(:,142)+ALGEBRAIC(:,46)+ALGEBRAIC(:,110)+ALGEBRAIC(:,75)+ALGEBRAIC(:,108)+ CONSTANTS(:,1).*ALGEBRAIC(:,118)+ALGEBRAIC(:,50)+ALGEBRAIC(:,83)) -  CONSTANTS(:,1).*ALGEBRAIC(:,122))+ALGEBRAIC(:,54)+ALGEBRAIC(:,128)+ALGEBRAIC(:,114)+ALGEBRAIC(:,73)+ 2.00000.*CONSTANTS(:,1).*ALGEBRAIC(:,120)+ALGEBRAIC(:,138)) - ALGEBRAIC(:,58)) - ALGEBRAIC(:,68);
    RATES(:,25) =  - ALGEBRAIC(:,145)./( CONSTANTS(:,75).*CONSTANTS(:,106));
   RATES = RATES';
end

% Calculate algebraic variables
function ALGEBRAIC = computeAlgebraic(ALGEBRAIC, CONSTANTS, STATES, VOI)
    statesSize = size(STATES);
    statesColumnCount = statesSize(2);
    if ( statesColumnCount == 1)
        STATES = STATES';
        utilOnes = 1;
    else
        statesRowCount = statesSize(1);
        utilOnes = ones(statesRowCount, 1);
    end
    ALGEBRAIC(:,1) = (1.00000+power(CONSTANTS(:,9)./STATES(:,7), 4.00000)+power(STATES(:,7)./CONSTANTS(:,10), 3.00000))./(1.00000+1.00000./CONSTANTS(:,11)+power(CONSTANTS(:,9)./STATES(:,7), 4.00000)+power(STATES(:,7)./CONSTANTS(:,10), 3.00000));
    ALGEBRAIC(:,2) = (1.00000+power(CONSTANTS(:,9)./STATES(:,15), 4.00000)+power(STATES(:,15)./CONSTANTS(:,10), 3.00000))./(1.00000+1.00000./CONSTANTS(:,11)+power(CONSTANTS(:,9)./STATES(:,15), 4.00000)+power(STATES(:,15)./CONSTANTS(:,10), 3.00000));
    ALGEBRAIC(:,9) = 1.00000./(1.00000+exp( - (STATES(:,24)+0.0239000)./0.00480000));
    ALGEBRAIC(:,10) = 1.00000./(1.00000+exp( - (STATES(:,25)+0.0239000)./0.00480000));
    ALGEBRAIC(:,11) = 1.00000./(1.00000+exp((STATES(:,24)+0.0661000)./0.00650000));
    ALGEBRAIC(:,12) = 1.00000./(1.00000+exp((STATES(:,25)+0.0661000)./0.00650000));
    ALGEBRAIC(:,18) = (  - CONSTANTS(:,1).*STATES(:,24))./( CONSTANTS(:,76).*CONSTANTS(:,77));
    ALGEBRAIC(:,20) =  CONSTANTS(:,97).*exp(CONSTANTS(:,95)+ CONSTANTS(:,96).*STATES(:,24));
    ALGEBRAIC(:,21) =  CONSTANTS(:,97).*exp(CONSTANTS(:,95)+ CONSTANTS(:,96).*STATES(:,25));
    ALGEBRAIC(:,5) = 1.00000./(1.00000+exp( - (STATES(:,24)+0.00177000)./0.0145200));
    ALGEBRAIC(:,22) =  0.210987.*exp( - power((STATES(:,24)+0.214340)./0.195350, 2.00000)) - 0.0205900;
    ALGEBRAIC(:,6) = 1.00000./(1.00000+exp( - (STATES(:,25)+0.00177000)./0.0145200));
    ALGEBRAIC(:,23) =  0.210987.*exp( - power((STATES(:,25)+0.214340)./0.195350, 2.00000)) - 0.0205900;
    ALGEBRAIC(:,24) =  0.821390.*exp( - power((STATES(:,24)+0.0315900)./0.0274600, 2.00000))+0.000190000;
    ALGEBRAIC(:,25) =  0.821390.*exp( - power((STATES(:,25)+0.0315900)./0.0274600, 2.00000))+0.000190000;
    ALGEBRAIC(:,7) =  ((( CONSTANTS(:,134).*arbitrary_log(STATES(:,7)./1.00000, 10) - CONSTANTS(:,135)) -  CONSTANTS(:,133).*(STATES(:,69)./(CONSTANTS(:,132)+STATES(:,69)))) -  CONSTANTS(:,131).*(power(STATES(:,67), 2.00000)./(power(CONSTANTS(:,130), 2.00000)+power(STATES(:,67), 2.00000)))).*CONSTANTS(:,171);
    ALGEBRAIC(:,26) = 1.00000./(1.00000+exp( - (STATES(:,24) - ALGEBRAIC(:,7))./0.0308000));
    ALGEBRAIC(:,8) =  ((( CONSTANTS(:,134).*arbitrary_log(STATES(:,15)./1.00000, 10) - CONSTANTS(:,135)) -  CONSTANTS(:,133).*(STATES(:,69)./(CONSTANTS(:,132)+STATES(:,69)))) -  CONSTANTS(:,131).*(power(STATES(:,67), 2.00000)./(power(CONSTANTS(:,130), 2.00000)+power(STATES(:,67), 2.00000)))).*CONSTANTS(:,171);
    ALGEBRAIC(:,27) = 1.00000./(1.00000+exp( - (STATES(:,25) - ALGEBRAIC(:,8))./0.0308000));
    ALGEBRAIC(:,13) = 1.00000./(1.00000+exp( - (STATES(:,24)+0.00188000)./0.00757000));
    ALGEBRAIC(:,28) =  0.00289000.*exp( - power((STATES(:,24)+0.00863000)./0.0123900, 2.00000))+0.00243000;
    ALGEBRAIC(:,14) = 1.00000./(1.00000+exp( - (STATES(:,25)+0.00188000)./0.00757000));
    ALGEBRAIC(:,29) =  0.0289000.*exp( - power((STATES(:,25)+0.00863000)./0.0123900, 2.00000))+0.00243000;
    ALGEBRAIC(:,15) = 1.00000./(1.00000+exp((STATES(:,24)+0.0293200)./0.00154000));
    ALGEBRAIC(:,30) =  0.295590.*exp( - power((STATES(:,24) - 0.00472000)./0.112550, 2.00000))+0.0231900;
    ALGEBRAIC(:,16) = 1.00000./(1.00000+exp((STATES(:,25)+0.0293200)./0.00154000));
    ALGEBRAIC(:,31) =  0.295590.*exp( - power((STATES(:,25) - 0.00472000)./0.112550, 2.00000))+0.0231900;
    ALGEBRAIC(:,19) = 1.00000 - (STATES(:,71)+STATES(:,72));
    ALGEBRAIC(:,33) =  CONSTANTS(:,123).*STATES(:,67);
    ALGEBRAIC(:,32) =  CONSTANTS(:,157).*STATES(:,12);
    ALGEBRAIC(:,36) = 1.00000+ CONSTANTS(:,142).*(power(STATES(:,68), 2.00000)./(power(CONSTANTS(:,141), 2.00000)+power(STATES(:,68), 2.00000)));
    ALGEBRAIC(:,39) = ( CONSTANTS(:,158).*(1.00000+ 3.65000.*(power(STATES(:,67), 2.00000)./(power(CONSTANTS(:,170), 2.00000)+power(STATES(:,67), 2.00000)))))./ALGEBRAIC(:,36);
    ALGEBRAIC(:,17) = CONSTANTS(:,164) - (STATES(:,64)+STATES(:,65)+STATES(:,66));
    ALGEBRAIC(:,41) = ALGEBRAIC(:,39);
    ALGEBRAIC(:,63) =  (STATES(:,28)./CONSTANTS(:,19)).*STATES(:,34);
    ALGEBRAIC(:,61) =  (STATES(:,28)./CONSTANTS(:,91)).*STATES(:,32);
    ALGEBRAIC(:,59) =  (STATES(:,28)./CONSTANTS(:,19)).*STATES(:,30);
    ALGEBRAIC(:,67) =  (STATES(:,29)./CONSTANTS(:,19)).*STATES(:,35);
    ALGEBRAIC(:,62) =  (STATES(:,29)./CONSTANTS(:,91)).*STATES(:,33);
    ALGEBRAIC(:,60) =  (STATES(:,29)./CONSTANTS(:,19)).*STATES(:,31);
    ALGEBRAIC(:,35) =  (CONSTANTS(:,179)./CONSTANTS(:,87)).*log(CONSTANTS(:,83)./STATES(:,1));
    ALGEBRAIC(:,47) =  (CONSTANTS(:,179)./CONSTANTS(:,89)).*log(CONSTANTS(:,85)./STATES(:,5));
    ALGEBRAIC(:,78) =  CONSTANTS(:,28).*CONSTANTS(:,74).*((ALGEBRAIC(:,35) - ALGEBRAIC(:,47))./((ALGEBRAIC(:,35) - ALGEBRAIC(:,47))+CONSTANTS(:,26)));
    ALGEBRAIC(:,80) =  CONSTANTS(:,29).*CONSTANTS(:,74).*(( ( CONSTANTS(:,84).*CONSTANTS(:,83).*power(CONSTANTS(:,85), 2.00000) -  STATES(:,3).*STATES(:,1).*power(STATES(:,5), 2.00000)).*CONSTANTS(:,30).*CONSTANTS(:,31).*CONSTANTS(:,32).*CONSTANTS(:,33))./( (CONSTANTS(:,84)+CONSTANTS(:,30)).*(CONSTANTS(:,83)+CONSTANTS(:,31)).*(CONSTANTS(:,85)+CONSTANTS(:,32)).*(CONSTANTS(:,85)+CONSTANTS(:,33)).*(STATES(:,3)+CONSTANTS(:,30)).*(STATES(:,1)+CONSTANTS(:,31)).*(STATES(:,5)+CONSTANTS(:,32)).*(STATES(:,5)+CONSTANTS(:,33))));
    ALGEBRAIC(:,43) =  (CONSTANTS(:,179)./CONSTANTS(:,88)).*log(CONSTANTS(:,84)./STATES(:,3));
    ALGEBRAIC(:,76) =  CONSTANTS(:,27).*CONSTANTS(:,74).*(power(ALGEBRAIC(:,43) - ALGEBRAIC(:,47), 4.00000)./(power(ALGEBRAIC(:,43) - ALGEBRAIC(:,47), 4.00000)+power(CONSTANTS(:,26), 4.00000)));
    ALGEBRAIC(:,48) =  CONSTANTS(:,43).*CONSTANTS(:,74).*(STATES(:,24) - ALGEBRAIC(:,47));
    ALGEBRAIC(:,82) =  CONSTANTS(:,39).*CONSTANTS(:,74).*(STATES(:,39)+STATES(:,40)+STATES(:,41)).*(STATES(:,24) - ALGEBRAIC(:,47));
    ALGEBRAIC(:,121) =  (( CONSTANTS(:,104).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77))).*(STATES(:,25) - STATES(:,24));
    ALGEBRAIC(:,122) = piecewise({STATES(:,25) - STATES(:,24)==0.00000,  CONSTANTS(:,186).*(STATES(:,6) - STATES(:,5)) },  CONSTANTS(:,186).*ALGEBRAIC(:,121).*((STATES(:,6) -  STATES(:,5).*exp( - ALGEBRAIC(:,121)))./(1.00000 - exp( - ALGEBRAIC(:,121)))));
    ALGEBRAIC(:,40) =  (CONSTANTS(:,179)./CONSTANTS(:,87)).*log(CONSTANTS(:,83)./STATES(:,2));
    ALGEBRAIC(:,49) =  (CONSTANTS(:,179)./CONSTANTS(:,89)).*log(CONSTANTS(:,85)./STATES(:,6));
    ALGEBRAIC(:,79) =  CONSTANTS(:,28).*CONSTANTS(:,75).*((ALGEBRAIC(:,40) - ALGEBRAIC(:,49))./((ALGEBRAIC(:,40) - ALGEBRAIC(:,49))+CONSTANTS(:,26)));
    ALGEBRAIC(:,81) =  CONSTANTS(:,29).*CONSTANTS(:,75).*(( ( CONSTANTS(:,84).*CONSTANTS(:,83).*power(CONSTANTS(:,85), 2.00000) -  STATES(:,4).*STATES(:,2).*power(STATES(:,6), 2.00000)).*CONSTANTS(:,30).*CONSTANTS(:,31).*CONSTANTS(:,32).*CONSTANTS(:,33))./( (CONSTANTS(:,84)+CONSTANTS(:,30)).*(CONSTANTS(:,83)+CONSTANTS(:,31)).*(CONSTANTS(:,85)+CONSTANTS(:,32)).*(CONSTANTS(:,85)+CONSTANTS(:,33)).*(STATES(:,4)+CONSTANTS(:,30)).*(STATES(:,2)+CONSTANTS(:,31)).*(STATES(:,6)+CONSTANTS(:,32)).*(STATES(:,6)+CONSTANTS(:,33))));
    ALGEBRAIC(:,45) =  (CONSTANTS(:,179)./CONSTANTS(:,88)).*log(CONSTANTS(:,84)./STATES(:,4));
    ALGEBRAIC(:,77) =  CONSTANTS(:,27).*CONSTANTS(:,75).*(power(ALGEBRAIC(:,45) - ALGEBRAIC(:,49), 4.00000)./(power(ALGEBRAIC(:,45) - ALGEBRAIC(:,49), 4.00000)+power(CONSTANTS(:,26), 4.00000)));
    ALGEBRAIC(:,50) =  CONSTANTS(:,43).*CONSTANTS(:,75).*(STATES(:,25) - ALGEBRAIC(:,49));
    ALGEBRAIC(:,83) =  CONSTANTS(:,39).*CONSTANTS(:,75).*(STATES(:,45)+STATES(:,46)+STATES(:,47)).*(STATES(:,25) - ALGEBRAIC(:,49));
    ALGEBRAIC(:,123) = CONSTANTS(:,162) - (STATES(:,13)+STATES(:,8)+STATES(:,10)+STATES(:,11)+STATES(:,12));
    ALGEBRAIC(:,124) = CONSTANTS(:,163) - (STATES(:,11)+STATES(:,12));
    ALGEBRAIC(:,129) = power(STATES(:,70), 1.0 ./ 2)./(power(CONSTANTS(:,137), 1.0 ./ 2)+power(STATES(:,70), 1.0 ./ 2));
    ALGEBRAIC(:,130) = STATES(:,68)./(CONSTANTS(:,138)+STATES(:,68));
    ALGEBRAIC(:,131) =  CONSTANTS(:,167).*(1.00000 - ALGEBRAIC(:,129)).*(1.00000 - ALGEBRAIC(:,130));
    ALGEBRAIC(:,125) =  (CONSTANTS(:,165)./CONSTANTS(:,172)).*(1.00000 - STATES(:,67)./( 2.00000.*(CONSTANTS(:,128)+STATES(:,67))));
    ALGEBRAIC(:,132) =  ALGEBRAIC(:,131).*CONSTANTS(:,74).*((power(STATES(:,7)./ALGEBRAIC(:,125), CONSTANTS(:,7)) - power(STATES(:,20)./CONSTANTS(:,6), CONSTANTS(:,7)))./(1.00000+(power(STATES(:,7)./ALGEBRAIC(:,125), CONSTANTS(:,7))+power(STATES(:,20)./CONSTANTS(:,6), CONSTANTS(:,7)))));
    ALGEBRAIC(:,64) =  CONSTANTS(:,13).*power(ALGEBRAIC(:,63), 3.00000).*(STATES(:,20) - STATES(:,7)).*2.00000.*CONSTANTS(:,1).*CONSTANTS(:,4);
    ALGEBRAIC(:,68) =  CONSTANTS(:,13).*power(ALGEBRAIC(:,67), 3.00000).*(STATES(:,20) - STATES(:,15)).*2.00000.*CONSTANTS(:,1).*CONSTANTS(:,4);
    ALGEBRAIC(:,133) =  ALGEBRAIC(:,131).*CONSTANTS(:,75).*((power(STATES(:,15)./ALGEBRAIC(:,125), CONSTANTS(:,7)) - power(STATES(:,20)./CONSTANTS(:,6), CONSTANTS(:,7)))./(1.00000+(power(STATES(:,15)./ALGEBRAIC(:,125), CONSTANTS(:,7))+power(STATES(:,20)./CONSTANTS(:,6), CONSTANTS(:,7)))));
    ALGEBRAIC(:,51) =  (CONSTANTS(:,179)./CONSTANTS(:,90)).*log(CONSTANTS(:,86)./STATES(:,7));
    ALGEBRAIC(:,52) =  CONSTANTS(:,42).*CONSTANTS(:,74).*(STATES(:,24) - ALGEBRAIC(:,51));
    ALGEBRAIC(:,126) =  CONSTANTS(:,166).*(1.00000+( 3.60000.*STATES(:,67))./(CONSTANTS(:,129)+STATES(:,67)));
    ALGEBRAIC(:,127) =  ALGEBRAIC(:,126).*CONSTANTS(:,74).*(STATES(:,7)./(CONSTANTS(:,66)+STATES(:,7)));
    ALGEBRAIC(:,71) = (STATES(:,20)+STATES(:,21))./2.00000;
    ALGEBRAIC(:,72) =  CONSTANTS(:,24).*(CONSTANTS(:,25)./(CONSTANTS(:,25)+ALGEBRAIC(:,71))).*(STATES(:,24) - ALGEBRAIC(:,51));
    ALGEBRAIC(:,111) =  0.740000.*STATES(:,62)+0.260000;
    ALGEBRAIC(:,112) =  CONSTANTS(:,67).*CONSTANTS(:,74).*STATES(:,60).*ALGEBRAIC(:,111).*(STATES(:,24) - ALGEBRAIC(:,51));
    ALGEBRAIC(:,55) = ( STATES(:,26).*(1.00000+power(STATES(:,7)./CONSTANTS(:,10), 3.00000)))./(1.00000+power(CONSTANTS(:,9)./STATES(:,7), 4.00000)+power(STATES(:,7)./CONSTANTS(:,10), 3.00000));
    ALGEBRAIC(:,56) =  CONSTANTS(:,8).*ALGEBRAIC(:,55).*(STATES(:,21) - STATES(:,7)).*( 2.00000.*CONSTANTS(:,1).*CONSTANTS(:,4));
    ALGEBRAIC(:,119) =  (( CONSTANTS(:,102).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77))).*(STATES(:,25) - STATES(:,24));
    ALGEBRAIC(:,120) = piecewise({STATES(:,25) - STATES(:,24)==0.00000,  CONSTANTS(:,185).*(STATES(:,15) - STATES(:,7)) },  CONSTANTS(:,185).*ALGEBRAIC(:,119).*((STATES(:,15) -  STATES(:,7).*exp( - ALGEBRAIC(:,119)))./(1.00000 - exp( - ALGEBRAIC(:,119)))));
    ALGEBRAIC(:,134) =  ALGEBRAIC(:,131).*CONSTANTS(:,74).*((power(STATES(:,7)./ALGEBRAIC(:,125), CONSTANTS(:,7)) - power(STATES(:,21)./CONSTANTS(:,6), CONSTANTS(:,7)))./(1.00000+(power(STATES(:,7)./ALGEBRAIC(:,125), CONSTANTS(:,7))+power(STATES(:,21)./CONSTANTS(:,6), CONSTANTS(:,7)))));
    ALGEBRAIC(:,135) =  ALGEBRAIC(:,131).*CONSTANTS(:,75).*((power(STATES(:,15)./ALGEBRAIC(:,125), CONSTANTS(:,7)) - power(STATES(:,21)./CONSTANTS(:,6), CONSTANTS(:,7)))./(1.00000+(power(STATES(:,15)./ALGEBRAIC(:,125), CONSTANTS(:,7))+power(STATES(:,21)./CONSTANTS(:,6), CONSTANTS(:,7)))));
    ALGEBRAIC(:,57) = ( STATES(:,27).*(1.00000+power(STATES(:,15)./CONSTANTS(:,10), 3.00000)))./(1.00000+power(CONSTANTS(:,9)./STATES(:,15), 4.00000)+power(STATES(:,15)./CONSTANTS(:,10), 3.00000));
    ALGEBRAIC(:,58) =  CONSTANTS(:,8).*ALGEBRAIC(:,57).*(STATES(:,21) - STATES(:,15)).*( 2.00000.*CONSTANTS(:,1).*CONSTANTS(:,4));
    ALGEBRAIC(:,137) = CONSTANTS(:,162) - (STATES(:,73)+STATES(:,16)+STATES(:,17)+STATES(:,18)+STATES(:,19));
    ALGEBRAIC(:,140) = CONSTANTS(:,163) - (STATES(:,18)+STATES(:,19));
    ALGEBRAIC(:,38) =  CONSTANTS(:,40).*CONSTANTS(:,74).*(STATES(:,24) - ALGEBRAIC(:,35));
    ALGEBRAIC(:,84) = 3.24100./(1.00000+exp(((STATES(:,24) - ALGEBRAIC(:,35)) - 185.720)./1.59600));
    ALGEBRAIC(:,86) =  1.00000.*(( 0.0130000.*exp(((STATES(:,24) - ALGEBRAIC(:,35))+5.79220)./0.294000)+exp(((STATES(:,24) - ALGEBRAIC(:,35)) - 653.733)./0.244000))./(1.00000+exp( - ((STATES(:,24) - ALGEBRAIC(:,35)) - 7.77750)./0.0928000)));
    ALGEBRAIC(:,87) = ALGEBRAIC(:,84)./(ALGEBRAIC(:,84)+ALGEBRAIC(:,86));
    ALGEBRAIC(:,88) =  CONSTANTS(:,44).*CONSTANTS(:,74).*ALGEBRAIC(:,87).*power((CONSTANTS(:,83)./CONSTANTS(:,45)), 1.0 ./ 2).*(STATES(:,24) - ALGEBRAIC(:,35));
    ALGEBRAIC(:,96) =  CONSTANTS(:,48).*CONSTANTS(:,74).*power(CONSTANTS(:,83)./CONSTANTS(:,45), CONSTANTS(:,47)).*(STATES(:,24) - ALGEBRAIC(:,35));
    ALGEBRAIC(:,92) =  0.580000.*STATES(:,48)+ 0.420000.*STATES(:,50);
    ALGEBRAIC(:,93) =  CONSTANTS(:,46).*CONSTANTS(:,74).*power(ALGEBRAIC(:,92), 2.00000).*(STATES(:,24) - ALGEBRAIC(:,35));
    ALGEBRAIC(:,98) =  0.650000.*STATES(:,52)+ 0.350000.*STATES(:,54);
    ALGEBRAIC(:,99) =  CONSTANTS(:,49).*CONSTANTS(:,74).*ALGEBRAIC(:,98).*(STATES(:,24) - ALGEBRAIC(:,35));
    ALGEBRAIC(:,139) =  CONSTANTS(:,168).*(1.00000 - STATES(:,70)./(CONSTANTS(:,139)+STATES(:,70))).*(1.00000 - STATES(:,68)./(CONSTANTS(:,140)+STATES(:,68)));
    ALGEBRAIC(:,102) = 1.00000./(1.00000+ 0.124500.*exp((  - 0.100000.*STATES(:,24).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77)))+ 0.0365000.*CONSTANTS(:,180).*exp((  - STATES(:,24).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77))));
    ALGEBRAIC(:,141) =  ALGEBRAIC(:,139).*ALGEBRAIC(:,102).*(CONSTANTS(:,83)./(CONSTANTS(:,83)+CONSTANTS(:,52))).*(power(STATES(:,3), 1.50000)./(power(STATES(:,3), 1.50000)+power(CONSTANTS(:,53), 1.50000)));
    ALGEBRAIC(:,115) =  (( CONSTANTS(:,98).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77))).*(STATES(:,25) - STATES(:,24));
    ALGEBRAIC(:,116) = piecewise({STATES(:,25) - STATES(:,24)==0.00000,  CONSTANTS(:,187).*(STATES(:,2) - STATES(:,1)) },  CONSTANTS(:,187).*ALGEBRAIC(:,115).*((STATES(:,2) -  STATES(:,1).*exp( - ALGEBRAIC(:,115)))./(1.00000 - exp( - ALGEBRAIC(:,115)))));
    ALGEBRAIC(:,42) =  CONSTANTS(:,40).*CONSTANTS(:,75).*(STATES(:,25) - ALGEBRAIC(:,40));
    ALGEBRAIC(:,85) = 3.24100./(1.00000+exp(((STATES(:,25) - ALGEBRAIC(:,40)) - 185.720)./1.59600));
    ALGEBRAIC(:,89) =  1.00000.*(( 0.0130000.*exp(((STATES(:,25) - ALGEBRAIC(:,40))+5.79220)./0.294000)+exp(((STATES(:,25) - ALGEBRAIC(:,40)) - 653.733)./0.244000))./(1.00000+exp( - (((STATES(:,25) - ALGEBRAIC(:,40)) - 7.77750)./0.0928000))));
    ALGEBRAIC(:,90) = ALGEBRAIC(:,85)./(ALGEBRAIC(:,85)+ALGEBRAIC(:,89));
    ALGEBRAIC(:,91) =  CONSTANTS(:,44).*CONSTANTS(:,75).*ALGEBRAIC(:,90).*power((CONSTANTS(:,83)./CONSTANTS(:,45)), 1.0 ./ 2).*(STATES(:,25) - ALGEBRAIC(:,40));
    ALGEBRAIC(:,97) =  CONSTANTS(:,48).*CONSTANTS(:,75).*power(CONSTANTS(:,83)./CONSTANTS(:,45), CONSTANTS(:,47)).*(STATES(:,25) - ALGEBRAIC(:,40));
    ALGEBRAIC(:,94) =  0.580000.*STATES(:,49)+ 0.420000.*STATES(:,51);
    ALGEBRAIC(:,95) =  CONSTANTS(:,46).*CONSTANTS(:,75).*power(ALGEBRAIC(:,94), 2.00000).*(STATES(:,25) - ALGEBRAIC(:,40));
    ALGEBRAIC(:,100) =  0.650000.*STATES(:,53)+ 0.350000.*STATES(:,55);
    ALGEBRAIC(:,101) =  CONSTANTS(:,49).*CONSTANTS(:,75).*ALGEBRAIC(:,100).*(STATES(:,25) - ALGEBRAIC(:,40));
    ALGEBRAIC(:,103) = 1.00000./(1.00000+ 0.124500.*exp((  - 0.100000.*STATES(:,25).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77)))+ 0.0365000.*CONSTANTS(:,180).*exp((  - STATES(:,25).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77))));
    ALGEBRAIC(:,142) =  ALGEBRAIC(:,139).*ALGEBRAIC(:,103).*(CONSTANTS(:,83)./(CONSTANTS(:,83)+CONSTANTS(:,52))).*(power(STATES(:,4), 1.50000)./(power(STATES(:,4), 1.50000)+power(CONSTANTS(:,54), 1.50000)));
    ALGEBRAIC(:,44) =  CONSTANTS(:,41).*CONSTANTS(:,74).*(STATES(:,24) - ALGEBRAIC(:,43));
    ALGEBRAIC(:,109) =  CONSTANTS(:,63).*CONSTANTS(:,74).*STATES(:,56).*STATES(:,58).*(STATES(:,24) - ALGEBRAIC(:,43));
    ALGEBRAIC(:,74) =  ALGEBRAIC(:,72).*( (power(CONSTANTS(:,92), 2.00000)./power(CONSTANTS(:,93), 2.00000)).*CONSTANTS(:,94)).*((STATES(:,3) -  CONSTANTS(:,84).*exp((  - CONSTANTS(:,92).*CONSTANTS(:,1).*STATES(:,24))./( CONSTANTS(:,76).*CONSTANTS(:,77))))./(STATES(:,7) -  CONSTANTS(:,86).*exp((  - CONSTANTS(:,93).*CONSTANTS(:,1).*STATES(:,24))./( CONSTANTS(:,76).*CONSTANTS(:,77))))).*((1.00000 - exp((  - CONSTANTS(:,93).*CONSTANTS(:,1).*STATES(:,24))./( CONSTANTS(:,76).*CONSTANTS(:,77))))./(1.00000 - exp((  - CONSTANTS(:,92).*CONSTANTS(:,1).*STATES(:,24))./( CONSTANTS(:,76).*CONSTANTS(:,77)))));
    ALGEBRAIC(:,117) =  (( CONSTANTS(:,100).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77))).*(STATES(:,25) - STATES(:,24));
    ALGEBRAIC(:,118) = piecewise({STATES(:,25) - STATES(:,24)==0.00000,  CONSTANTS(:,188).*(STATES(:,4) - STATES(:,3)) },  CONSTANTS(:,188).*ALGEBRAIC(:,117).*((STATES(:,4) -  STATES(:,3).*exp( - ALGEBRAIC(:,117)))./(1.00000 - exp( - ALGEBRAIC(:,117)))));
    ALGEBRAIC(:,46) =  CONSTANTS(:,41).*CONSTANTS(:,75).*(STATES(:,25) - ALGEBRAIC(:,45));
    ALGEBRAIC(:,110) =  CONSTANTS(:,63).*CONSTANTS(:,75).*STATES(:,57).*STATES(:,59).*(STATES(:,25) - ALGEBRAIC(:,45));
    ALGEBRAIC(:,53) =  (CONSTANTS(:,179)./CONSTANTS(:,90)).*log(CONSTANTS(:,86)./STATES(:,15));
    ALGEBRAIC(:,73) =  CONSTANTS(:,24).*(CONSTANTS(:,25)./(CONSTANTS(:,25)+ALGEBRAIC(:,71))).*(STATES(:,25) - ALGEBRAIC(:,53));
    ALGEBRAIC(:,75) =  ALGEBRAIC(:,73).*( (power(CONSTANTS(:,92), 2.00000)./power(CONSTANTS(:,93), 2.00000)).*CONSTANTS(:,94)).*((STATES(:,4) -  CONSTANTS(:,84).*exp((  - CONSTANTS(:,92).*CONSTANTS(:,1).*STATES(:,25))./( CONSTANTS(:,76).*CONSTANTS(:,77))))./(STATES(:,15) -  CONSTANTS(:,86).*exp((  - CONSTANTS(:,93).*CONSTANTS(:,1).*STATES(:,25))./( CONSTANTS(:,76).*CONSTANTS(:,77))))).*((1.00000 - exp((  - CONSTANTS(:,93).*CONSTANTS(:,1).*STATES(:,25))./( CONSTANTS(:,76).*CONSTANTS(:,77))))./(1.00000 - exp((  - CONSTANTS(:,92).*CONSTANTS(:,1).*STATES(:,25))./( CONSTANTS(:,76).*CONSTANTS(:,77)))));
    ALGEBRAIC(:,104) = 1.00000./(1.00000+power(CONSTANTS(:,62)./STATES(:,15), 2.00000));
    ALGEBRAIC(:,105) = exp(( CONSTANTS(:,61).*STATES(:,25).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77)));
    ALGEBRAIC(:,106) = exp(( (CONSTANTS(:,61) - 1.00000).*STATES(:,25).*CONSTANTS(:,1))./( CONSTANTS(:,76).*CONSTANTS(:,77)));
    ALGEBRAIC(:,107) =  power(CONSTANTS(:,84), 3.00000).*STATES(:,15)+ power(STATES(:,4), 3.00000).*CONSTANTS(:,86)+ power(CONSTANTS(:,59), 3.00000).*STATES(:,15)+ CONSTANTS(:,57).*power(STATES(:,4), 3.00000)+ power(CONSTANTS(:,58), 3.00000).*CONSTANTS(:,86).*(1.00000+STATES(:,15)./CONSTANTS(:,56))+ CONSTANTS(:,56).*power(CONSTANTS(:,84), 3.00000).*(1.00000+power(STATES(:,4)./CONSTANTS(:,58), 3.00000));
    ALGEBRAIC(:,108) =  CONSTANTS(:,55).*ALGEBRAIC(:,104).*(( ALGEBRAIC(:,105).*power(STATES(:,4), 3.00000).*CONSTANTS(:,86) -  ALGEBRAIC(:,106).*power(CONSTANTS(:,84), 3.00000).*STATES(:,15))./( ALGEBRAIC(:,107).*(1.00000+ CONSTANTS(:,60).*ALGEBRAIC(:,106))));
    ALGEBRAIC(:,54) =  CONSTANTS(:,42).*CONSTANTS(:,75).*(STATES(:,25) - ALGEBRAIC(:,53));
    ALGEBRAIC(:,128) =  ALGEBRAIC(:,126).*CONSTANTS(:,75).*(STATES(:,15)./(CONSTANTS(:,66)+STATES(:,15)));
    ALGEBRAIC(:,113) =  0.740000.*STATES(:,63)+0.260000;
    ALGEBRAIC(:,114) =  CONSTANTS(:,67).*CONSTANTS(:,75).*STATES(:,61).*ALGEBRAIC(:,113).*(STATES(:,25) - ALGEBRAIC(:,53));
    ALGEBRAIC(:,136) = ALGEBRAIC(:,132)+ALGEBRAIC(:,134);
    ALGEBRAIC(:,144) = ((((((((ALGEBRAIC(:,38)+ALGEBRAIC(:,88)+ALGEBRAIC(:,96)+ALGEBRAIC(:,93)+ALGEBRAIC(:,99)) -  CONSTANTS(:,1).*ALGEBRAIC(:,116))+ALGEBRAIC(:,141)+ALGEBRAIC(:,44)+ALGEBRAIC(:,109)+ALGEBRAIC(:,74)) -  CONSTANTS(:,1).*ALGEBRAIC(:,118))+ALGEBRAIC(:,48)+ALGEBRAIC(:,82)+ CONSTANTS(:,1).*ALGEBRAIC(:,122)+ALGEBRAIC(:,52)+ALGEBRAIC(:,127)+ALGEBRAIC(:,112)+ALGEBRAIC(:,72)) -  2.00000.*CONSTANTS(:,1).*ALGEBRAIC(:,120))+ALGEBRAIC(:,136)) - ALGEBRAIC(:,56)) - ALGEBRAIC(:,64);
    ALGEBRAIC(:,138) = ALGEBRAIC(:,133)+ALGEBRAIC(:,135);
    ALGEBRAIC(:,145) = ((((ALGEBRAIC(:,42)+ALGEBRAIC(:,91)+ALGEBRAIC(:,97)+ALGEBRAIC(:,95)+ALGEBRAIC(:,101)+ CONSTANTS(:,1).*ALGEBRAIC(:,116)+ALGEBRAIC(:,142)+ALGEBRAIC(:,46)+ALGEBRAIC(:,110)+ALGEBRAIC(:,75)+ALGEBRAIC(:,108)+ CONSTANTS(:,1).*ALGEBRAIC(:,118)+ALGEBRAIC(:,50)+ALGEBRAIC(:,83)) -  CONSTANTS(:,1).*ALGEBRAIC(:,122))+ALGEBRAIC(:,54)+ALGEBRAIC(:,128)+ALGEBRAIC(:,114)+ALGEBRAIC(:,73)+ 2.00000.*CONSTANTS(:,1).*ALGEBRAIC(:,120)+ALGEBRAIC(:,138)) - ALGEBRAIC(:,58)) - ALGEBRAIC(:,68);
    ALGEBRAIC(:,34) = (STATES(:,65)+ALGEBRAIC(:,17))./CONSTANTS(:,164);
    ALGEBRAIC(:,37) = ((CONSTANTS(:,173) -  CONSTANTS(:,174).*ALGEBRAIC(:,34))+ CONSTANTS(:,175).*power(ALGEBRAIC(:,34), 2.00000)) -  CONSTANTS(:,176).*power(ALGEBRAIC(:,34), 3.00000);
    ALGEBRAIC(:,143) = ALGEBRAIC(:,139);
end

% Functions required for solving differential algebraic equation
function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_0(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN)
    CONSTANTS = CONSTANTS_IN;
    STATES = STATES_IN;
    ALGEBRAIC = ALGEBRAIC_IN;
    global initialGuess_0;
    if (length(initialGuess_0) ~= 1), initialGuess_0 = 0.1;, end
    options = optimset('Display', 'off', 'TolX', 1E-6);
    if length(VOI) == 1
        residualfn = @(algebraicCandidate)residualSN_0(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES);
        ALGEBRAIC(:,3) = fsolve(residualfn, initialGuess_0, options);
        initialGuess_0 = ALGEBRAIC(:,3);
    else
        SET_ALGEBRAIC(:,3) = logical(1);
        for i=1:length(VOI)
            residualfn = @(algebraicCandidate)residualSN_0(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:));
            TEMP_ALGEBRAIC(:,3) = fsolve(residualfn, initialGuess_0, options);
            ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC);
            initialGuess_0 = TEMP_ALGEBRAIC(:,3);
        end
    end
end

function resid = residualSN_0(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES)
    ALGEBRAIC(:,3) = algebraicCandidate;
    resid = (1.00000) - (ALGEBRAIC(:,3)+STATES(:,36)+STATES(:,37)+STATES(:,38)+STATES(:,39)+STATES(:,40)+STATES(:,41));
end

% Functions required for solving differential algebraic equation
function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_1(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN)
    CONSTANTS = CONSTANTS_IN;
    STATES = STATES_IN;
    ALGEBRAIC = ALGEBRAIC_IN;
    global initialGuess_1;
    if (length(initialGuess_1) ~= 1), initialGuess_1 = 0.1;, end
    options = optimset('Display', 'off', 'TolX', 1E-6);
    if length(VOI) == 1
        residualfn = @(algebraicCandidate)residualSN_1(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES);
        ALGEBRAIC(:,4) = fsolve(residualfn, initialGuess_1, options);
        initialGuess_1 = ALGEBRAIC(:,4);
    else
        SET_ALGEBRAIC(:,4) = logical(1);
        for i=1:length(VOI)
            residualfn = @(algebraicCandidate)residualSN_1(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:));
            TEMP_ALGEBRAIC(:,4) = fsolve(residualfn, initialGuess_1, options);
            ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC);
            initialGuess_1 = TEMP_ALGEBRAIC(:,4);
        end
    end
end

function resid = residualSN_1(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES)
    ALGEBRAIC(:,4) = algebraicCandidate;
    resid = (1.00000) - (ALGEBRAIC(:,4)+STATES(:,42)+STATES(:,43)+STATES(:,44)+STATES(:,45)+STATES(:,46)+STATES(:,47));
end

% Functions required for solving differential algebraic equation
function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_2(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN)
    ALGEBRAIC = ALGEBRAIC_IN;
    CONSTANTS = CONSTANTS_IN;
    STATES = STATES_IN;
    global initialGuess_2;
    if (length(initialGuess_2) ~= 2), initialGuess_2 = [0.1,0.1];, end
    options = optimset('Display', 'off', 'TolX', 1E-6);
    if length(VOI) == 1
        residualfn = @(algebraicCandidate)residualSN_2(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES);
        soln = fsolve(residualfn, initialGuess_2, options);
        initialGuess_2 = soln;
        ALGEBRAIC(:,65) = soln(1);
        ALGEBRAIC(:,66) = soln(2);
    else
        SET_ALGEBRAIC(:,65) = logical(1);
        SET_ALGEBRAIC(:,66) = logical(1);
        for i=1:length(VOI)
            residualfn = @(algebraicCandidate)residualSN_2(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:));
            soln = fsolve(residualfn, initialGuess_2, options);
            initialGuess_2 = soln;
            TEMP_ALGEBRAIC(:,65) = soln(1);
            TEMP_ALGEBRAIC(:,66) = soln(2);
            ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC);
        end
    end
end

function resid = residualSN_2(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES)
    ALGEBRAIC(:,65) = algebraicCandidate(1);
    ALGEBRAIC(:,66) = algebraicCandidate(2);
    resid(1) = ALGEBRAIC(:,66) -  (STATES(:,28)./CONSTANTS(:,91)).*ALGEBRAIC(:,65);
    resid(2) = ALGEBRAIC(:,65) - (1.00000 - (STATES(:,30)+STATES(:,32)+STATES(:,34)+ALGEBRAIC(:,59)+ALGEBRAIC(:,61)+ALGEBRAIC(:,63)+ALGEBRAIC(:,66)));
end

% Functions required for solving differential algebraic equation
function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_3(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN)
    ALGEBRAIC = ALGEBRAIC_IN;
    CONSTANTS = CONSTANTS_IN;
    STATES = STATES_IN;
    global initialGuess_3;
    if (length(initialGuess_3) ~= 2), initialGuess_3 = [0.1,0.1];, end
    options = optimset('Display', 'off', 'TolX', 1E-6);
    if length(VOI) == 1
        residualfn = @(algebraicCandidate)residualSN_3(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES);
        soln = fsolve(residualfn, initialGuess_3, options);
        initialGuess_3 = soln;
        ALGEBRAIC(:,69) = soln(1);
        ALGEBRAIC(:,70) = soln(2);
    else
        SET_ALGEBRAIC(:,69) = logical(1);
        SET_ALGEBRAIC(:,70) = logical(1);
        for i=1:length(VOI)
            residualfn = @(algebraicCandidate)residualSN_3(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:));
            soln = fsolve(residualfn, initialGuess_3, options);
            initialGuess_3 = soln;
            TEMP_ALGEBRAIC(:,69) = soln(1);
            TEMP_ALGEBRAIC(:,70) = soln(2);
            ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC);
        end
    end
end

function resid = residualSN_3(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES)
    ALGEBRAIC(:,69) = algebraicCandidate(1);
    ALGEBRAIC(:,70) = algebraicCandidate(2);
    resid(1) = ALGEBRAIC(:,70) -  (STATES(:,29)./CONSTANTS(:,91)).*ALGEBRAIC(:,69);
    resid(2) = ALGEBRAIC(:,69) - (1.00000 - (STATES(:,31)+STATES(:,33)+STATES(:,35)+ALGEBRAIC(:,60)+ALGEBRAIC(:,62)+ALGEBRAIC(:,67)+ALGEBRAIC(:,70)));
end

% Compute result of a piecewise function
function x = piecewise(cases, default)
    set = [0];
    for i = 1:2:length(cases)
        if (length(cases{i+1}) == 1)
            x(cases{i} & ~set,:) = cases{i+1};
        else
            x(cases{i} & ~set,:) = cases{i+1}(cases{i} & ~set);
        end
        set = set | cases{i};
        if(set), break, end
    end
    if (length(default) == 1)
        x(~set,:) = default;
    else
        x(~set,:) = default(~set);
    end
end

% Compute a logarithm to any base" +
function x = arbitrary_log(a, base)
    x = log(a) ./ log(base);
end

% Pad out or shorten strings to a set length
function strout = strpad(strin)
    req_length = 160;
    insize = size(strin,2);
    if insize > req_length
        strout = strin(1:req_length);
    else
        strout = [strin, blanks(req_length - insize)];
    end
end