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 =126;
end
% There are a total of 41 entries in each of the rate and state variable arrays.
% There are a total of 126 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_CONSTANTS(:,1) = strpad('Bmax_CaM in component buffca (mM)');
    LEGEND_CONSTANTS(:,119) = strpad('Bmax_SLhighjn in component buffca (mM)');
    LEGEND_CONSTANTS(:,123) = strpad('Bmax_SLhighsl in component buffca (mM)');
    LEGEND_CONSTANTS(:,120) = strpad('Bmax_SLlowjn in component buffca (mM)');
    LEGEND_CONSTANTS(:,124) = strpad('Bmax_SLlowsl in component buffca (mM)');
    LEGEND_CONSTANTS(:,88) = strpad('Bmax_SR in component buffca (mM)');
    LEGEND_CONSTANTS(:,2) = strpad('Bmax_TnChigh in component buffca (mM)');
    LEGEND_CONSTANTS(:,3) = strpad('Bmax_TnClow in component buffca (mM)');
    LEGEND_CONSTANTS(:,4) = strpad('Bmax_myosin in component buffca (mM)');
    LEGEND_STATES(:,1) = strpad('CaM in component buffca (mM)');
    LEGEND_STATES(:,2) = strpad('Ca_i in component calcium (mM)');
    LEGEND_STATES(:,3) = strpad('Ca_jn in component calcium (mM)');
    LEGEND_STATES(:,4) = strpad('Ca_sl in component calcium (mM)');
    LEGEND_ALGEBRAIC(:,38) = strpad('J_CaB_cytosol in component buffca (mM_per_ms)');
    LEGEND_ALGEBRAIC(:,36) = strpad('J_CaB_jn in component buffca (mM_per_ms)');
    LEGEND_ALGEBRAIC(:,37) = strpad('J_CaB_sl in component buffca (mM_per_ms)');
    LEGEND_CONSTANTS(:,5) = strpad('Mg_i in component ion (mM)');
    LEGEND_STATES(:,5) = strpad('Myoc in component buffca (mM)');
    LEGEND_STATES(:,6) = strpad('Myom in component buffca (mM)');
    LEGEND_STATES(:,7) = strpad('SLH_jn in component buffca (mM)');
    LEGEND_STATES(:,8) = strpad('SLH_sl in component buffca (mM)');
    LEGEND_STATES(:,9) = strpad('SLL_jn in component buffca (mM)');
    LEGEND_STATES(:,10) = strpad('SLL_sl in component buffca (mM)');
    LEGEND_STATES(:,11) = strpad('SRB in component buffca (mM)');
    LEGEND_STATES(:,12) = strpad('TnCHc in component buffca (mM)');
    LEGEND_STATES(:,13) = strpad('TnCHm in component buffca (mM)');
    LEGEND_STATES(:,14) = strpad('TnCL in component buffca (mM)');
    LEGEND_CONSTANTS(:,114) = strpad('Vjn in component geom (litre)');
    LEGEND_CONSTANTS(:,117) = strpad('Vmyo in component geom (litre)');
    LEGEND_CONSTANTS(:,121) = strpad('Vsl in component geom (litre)');
    LEGEND_CONSTANTS(:,6) = strpad('koff_cam in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,7) = strpad('koff_myoca in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,8) = strpad('koff_myomg in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,9) = strpad('koff_slh in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,10) = strpad('koff_sll in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,11) = strpad('koff_sr in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,12) = strpad('koff_tnchca in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,13) = strpad('koff_tnchmg in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,14) = strpad('koff_tncl in component buffca (mS_per_uF)');
    LEGEND_CONSTANTS(:,15) = strpad('kon_cam in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,16) = strpad('kon_myoca in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,17) = strpad('kon_myomg in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,18) = strpad('kon_slh in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,19) = strpad('kon_sll in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,20) = strpad('kon_sr in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,21) = strpad('kon_tnchca in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,22) = strpad('kon_tnchmg in component buffca (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,23) = strpad('kon_tncl in component buffca (per_mM_per_ms)');
    LEGEND_VOI = strpad('time in component engine (ms)');
    LEGEND_CONSTANTS(:,24) = strpad('Bmax_Na_jn in component buffna (mM)');
    LEGEND_CONSTANTS(:,25) = strpad('Bmax_Na_sl in component buffna (mM)');
    LEGEND_STATES(:,15) = strpad('NaB_jn in component buffna (mM)');
    LEGEND_STATES(:,16) = strpad('NaB_sl in component buffna (mM)');
    LEGEND_STATES(:,17) = strpad('Na_jn in component sodium (mM)');
    LEGEND_STATES(:,18) = strpad('Na_sl in component sodium (mM)');
    LEGEND_CONSTANTS(:,26) = strpad('koff_na in component buffna (mS_per_uF)');
    LEGEND_CONSTANTS(:,27) = strpad('kon_na in component buffna (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,125) = strpad('Bmax_Csqn in component calcium (mM)');
    LEGEND_CONSTANTS(:,28) = strpad('C in component cell (farad)');
    LEGEND_STATES(:,19) = strpad('Ca_sr in component calcium (mM)');
    LEGEND_STATES(:,20) = strpad('Csqn in component calcium (mM)');
    LEGEND_CONSTANTS(:,29) = strpad('F in component phys (C_per_mol)');
    LEGEND_ALGEBRAIC(:,88) = strpad('ICaB_jn in component icab (A_per_F)');
    LEGEND_ALGEBRAIC(:,93) = strpad('ICaB_sl in component icab (A_per_F)');
    LEGEND_ALGEBRAIC(:,40) = strpad('ICaL_Ca_jn in component ical (A_per_F)');
    LEGEND_ALGEBRAIC(:,42) = strpad('ICaL_Ca_sl in component ical (A_per_F)');
    LEGEND_ALGEBRAIC(:,89) = strpad('ICa_tot_jn in component calcium (A_per_F)');
    LEGEND_ALGEBRAIC(:,96) = strpad('ICa_tot_sl in component calcium (A_per_F)');
    LEGEND_ALGEBRAIC(:,71) = strpad('INaCa_jn in component inaca (A_per_F)');
    LEGEND_ALGEBRAIC(:,75) = strpad('INaCa_sl in component inaca (A_per_F)');
    LEGEND_ALGEBRAIC(:,82) = strpad('IpCa_jn in component ipca (A_per_F)');
    LEGEND_ALGEBRAIC(:,84) = strpad('IpCa_sl in component ipca (A_per_F)');
    LEGEND_CONSTANTS(:,91) = strpad('JCa_jnsl in component calcium (m3_per_s)');
    LEGEND_CONSTANTS(:,92) = strpad('JCa_slmyo in component calcium (m3_per_s)');
    LEGEND_ALGEBRAIC(:,90) = strpad('J_SRCarel in component ryr (mM_per_ms)');
    LEGEND_ALGEBRAIC(:,92) = strpad('J_SRleak in component ryr (mM_per_ms)');
    LEGEND_ALGEBRAIC(:,94) = strpad('J_serca in component ryr (mM_per_ms)');
    LEGEND_CONSTANTS(:,122) = strpad('Vsr in component geom (litre)');
    LEGEND_CONSTANTS(:,30) = strpad('koff_csqn in component calcium (mS_per_uF)');
    LEGEND_CONSTANTS(:,31) = strpad('kon_csqn in component calcium (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,32) = strpad('AF in component cell (dimensionless)');
    LEGEND_CONSTANTS(:,89) = strpad('Vcell in component geom (litre)');
    LEGEND_CONSTANTS(:,33) = strpad('cell_length in component geom (um)');
    LEGEND_CONSTANTS(:,34) = strpad('cell_radius in component geom (um)');
    LEGEND_CONSTANTS(:,35) = strpad('pi in component geom (dimensionless)');
    LEGEND_ALGEBRAIC(:,87) = strpad('ECa_jn in component nernst (mV)');
    LEGEND_ALGEBRAIC(:,91) = strpad('ECa_sl in component nernst (mV)');
    LEGEND_CONSTANTS(:,36) = strpad('Fjn in component junc (dimensionless)');
    LEGEND_CONSTANTS(:,115) = strpad('Fsl in component junc (dimensionless)');
    LEGEND_ALGEBRAIC(:,95) = strpad('ICaB in component icab (A_per_F)');
    LEGEND_STATES(:,21) = strpad('V in component membrane (mV)');
    LEGEND_CONSTANTS(:,37) = strpad('gCaB in component icab (mS_per_uF)');
    LEGEND_CONSTANTS(:,38) = strpad('Ca_o in component ion (mM)');
    LEGEND_CONSTANTS(:,93) = strpad('FRT in component phys (per_mV)');
    LEGEND_CONSTANTS(:,39) = strpad('Fjn_CaL in component junc (dimensionless)');
    LEGEND_CONSTANTS(:,94) = strpad('Fsl_CaL in component junc (dimensionless)');
    LEGEND_ALGEBRAIC(:,53) = strpad('ICaL in component ical (A_per_F)');
    LEGEND_ALGEBRAIC(:,43) = strpad('ICaL_Ca in component ical (A_per_F)');
    LEGEND_ALGEBRAIC(:,45) = strpad('ICaL_K in component ical (A_per_F)');
    LEGEND_ALGEBRAIC(:,51) = strpad('ICaL_Na in component ical (A_per_F)');
    LEGEND_ALGEBRAIC(:,47) = strpad('ICaL_Na_jn in component ical (A_per_F)');
    LEGEND_ALGEBRAIC(:,49) = strpad('ICaL_Na_sl in component ical (A_per_F)');
    LEGEND_STATES(:,22) = strpad('K_i in component potassium (mM)');
    LEGEND_CONSTANTS(:,40) = strpad('K_o in component ion (mM)');
    LEGEND_CONSTANTS(:,41) = strpad('Na_o in component ion (mM)');
    LEGEND_CONSTANTS(:,118) = strpad('Q in component phys (dimensionless)');
    LEGEND_CONSTANTS(:,42) = strpad('Q10CaL in component ical (dimensionless)');
    LEGEND_STATES(:,23) = strpad('d in component ical (mS_per_uF)');
    LEGEND_STATES(:,24) = strpad('f in component ical (mS_per_uF)');
    LEGEND_STATES(:,25) = strpad('fCaB_jn in component ical (mS_per_uF)');
    LEGEND_STATES(:,26) = strpad('fCaB_sl in component ical (mS_per_uF)');
    LEGEND_CONSTANTS(:,43) = strpad('fcaCaMSL in component ical (dimensionless)');
    LEGEND_CONSTANTS(:,44) = strpad('fcaCaj in component ical (dimensionless)');
    LEGEND_ALGEBRAIC(:,39) = strpad('ibarca_jn in component ical (uA_per_cm2)');
    LEGEND_ALGEBRAIC(:,41) = strpad('ibarca_sl in component ical (uA_per_cm2)');
    LEGEND_ALGEBRAIC(:,44) = strpad('ibark in component ical (uA_per_cm2)');
    LEGEND_ALGEBRAIC(:,46) = strpad('ibarna_jn in component ical (uA_per_cm2)');
    LEGEND_ALGEBRAIC(:,48) = strpad('ibarna_sl in component ical (uA_per_cm2)');
    LEGEND_ALGEBRAIC(:,1) = strpad('ical_d_inf in component ical (dimensionless)');
    LEGEND_ALGEBRAIC(:,15) = strpad('ical_d_tau in component ical (ms)');
    LEGEND_ALGEBRAIC(:,2) = strpad('ical_f_inf in component ical (dimensionless)');
    LEGEND_ALGEBRAIC(:,16) = strpad('ical_f_tau in component ical (ms)');
    LEGEND_CONSTANTS(:,95) = strpad('pCa in component ical (cm_per_s)');
    LEGEND_CONSTANTS(:,96) = strpad('pK in component ical (cm_per_s)');
    LEGEND_CONSTANTS(:,97) = strpad('pNa in component ical (cm_per_s)');
    LEGEND_CONSTANTS(:,116) = strpad('ECl in component nernst (mV)');
    LEGEND_ALGEBRAIC(:,50) = strpad('IClB in component iclb (A_per_F)');
    LEGEND_CONSTANTS(:,45) = strpad('gClB in component iclb (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,55) = strpad('IClCa in component iclca (A_per_F)');
    LEGEND_ALGEBRAIC(:,52) = strpad('IClCa_jn in component iclca (A_per_F)');
    LEGEND_ALGEBRAIC(:,54) = strpad('IClCa_sl in component iclca (A_per_F)');
    LEGEND_CONSTANTS(:,46) = strpad('KdClCa in component iclca (mM)');
    LEGEND_CONSTANTS(:,47) = strpad('gClCa in component iclca (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,98) = strpad('EK in component nernst (mV)');
    LEGEND_ALGEBRAIC(:,102) = strpad('IK1 in component ik1 (A_per_F)');
    LEGEND_CONSTANTS(:,48) = strpad('IK1_Fact in component ik1 (dimensionless)');
    LEGEND_ALGEBRAIC(:,56) = strpad('Na_observed in component ik1 (mM)');
    LEGEND_ALGEBRAIC(:,57) = strpad('a_fact in component ik1 (dimensionless)');
    LEGEND_CONSTANTS(:,98) = strpad('gK1 in component ik1 (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,99) = strpad('ik1_inf_a in component ik1 (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,100) = strpad('ik1_inf_b in component ik1 (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,101) = strpad('inf in component ik1 (dimensionless)');
    LEGEND_CONSTANTS(:,99) = strpad('Act in component ikach (dimensionless)');
    LEGEND_CONSTANTS(:,49) = strpad('CCh in component ikach (dimensionless)');
    LEGEND_CONSTANTS(:,50) = strpad('Constitutive_Frac in component ikach (dimensionless)');
    LEGEND_ALGEBRAIC(:,103) = strpad('IKACh in component ikach (A_per_F)');
    LEGEND_ALGEBRAIC(:,59) = strpad('Na_dep in component ikach (dimensionless)');
    LEGEND_ALGEBRAIC(:,58) = strpad('Na_observed in component ikach (mM)');
    LEGEND_CONSTANTS(:,100) = strpad('gKACh in component ikach (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,106) = strpad('IKp in component ikp (A_per_F)');
    LEGEND_ALGEBRAIC(:,104) = strpad('IKp_jn in component ikp (A_per_F)');
    LEGEND_ALGEBRAIC(:,105) = strpad('IKp_sl in component ikp (A_per_F)');
    LEGEND_CONSTANTS(:,51) = strpad('gKp in component ikp (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,60) = strpad('kp in component ikp (dimensionless)');
    LEGEND_ALGEBRAIC(:,107) = strpad('IKr in component ikr (A_per_F)');
    LEGEND_CONSTANTS(:,101) = strpad('gKr in component ikr (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,3) = strpad('ikr_x_inf in component ikr (dimensionless)');
    LEGEND_ALGEBRAIC(:,17) = strpad('ikr_x_tau in component ikr (ms)');
    LEGEND_ALGEBRAIC(:,61) = strpad('r in component ikr (dimensionless)');
    LEGEND_STATES(:,27) = strpad('x in component ikr (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,62) = strpad('EKs in component iks (mV)');
    LEGEND_ALGEBRAIC(:,65) = strpad('IKs in component iks (A_per_F)');
    LEGEND_ALGEBRAIC(:,63) = strpad('IKs_jn in component iks (A_per_F)');
    LEGEND_ALGEBRAIC(:,64) = strpad('IKs_sl in component iks (A_per_F)');
    LEGEND_STATES(:,28) = strpad('Na_i in component sodium (mM)');
    LEGEND_CONSTANTS(:,102) = strpad('gKs_jn in component iks (mS_per_uF)');
    LEGEND_CONSTANTS(:,103) = strpad('gKs_sl in component iks (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,4) = strpad('iks_x_inf in component iks (dimensionless)');
    LEGEND_ALGEBRAIC(:,18) = strpad('iks_x_tau in component iks (ms)');
    LEGEND_CONSTANTS(:,52) = strpad('pNaK in component iks (dimensionless)');
    LEGEND_STATES(:,29) = strpad('x in component iks (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,108) = strpad('IKur in component ikur (A_per_F)');
    LEGEND_CONSTANTS(:,104) = strpad('gKur in component ikur (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,5) = strpad('ikur_r_inf in component ikur (dimensionless)');
    LEGEND_ALGEBRAIC(:,19) = strpad('ikur_r_tau in component ikur (ms)');
    LEGEND_ALGEBRAIC(:,6) = strpad('ikur_s_inf in component ikur (dimensionless)');
    LEGEND_ALGEBRAIC(:,20) = strpad('ikur_s_tau in component ikur (ms)');
    LEGEND_STATES(:,30) = strpad('r in component ikur (mS_per_uF)');
    LEGEND_STATES(:,31) = strpad('s in component ikur (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,110) = strpad('ENa_jn in component nernst (mV)');
    LEGEND_ALGEBRAIC(:,114) = strpad('ENa_sl in component nernst (mV)');
    LEGEND_ALGEBRAIC(:,116) = strpad('INa in component ina (A_per_F)');
    LEGEND_ALGEBRAIC(:,111) = strpad('INa_jn in component ina (A_per_F)');
    LEGEND_ALGEBRAIC(:,115) = strpad('INa_sl in component ina (A_per_F)');
    LEGEND_CONSTANTS(:,53) = strpad('gNa in component ina (mS_per_uF)');
    LEGEND_STATES(:,32) = strpad('h in component ina (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,7) = strpad('ina_h_alpha in component ina (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,21) = strpad('ina_h_beta in component ina (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,28) = strpad('ina_h_inf in component ina (dimensionless)');
    LEGEND_ALGEBRAIC(:,32) = strpad('ina_h_tau in component ina (ms)');
    LEGEND_ALGEBRAIC(:,8) = strpad('ina_j_alpha in component ina (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,22) = strpad('ina_j_beta in component ina (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,29) = strpad('ina_j_inf in component ina (dimensionless)');
    LEGEND_ALGEBRAIC(:,33) = strpad('ina_j_tau in component ina (ms)');
    LEGEND_ALGEBRAIC(:,9) = strpad('ina_m_alpha in component ina (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,23) = strpad('ina_m_beta in component ina (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,30) = strpad('ina_m_inf in component ina (dimensionless)');
    LEGEND_ALGEBRAIC(:,34) = strpad('ina_m_tau in component ina (ms)');
    LEGEND_STATES(:,33) = strpad('j in component ina (mS_per_uF)');
    LEGEND_STATES(:,34) = strpad('m in component ina (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,118) = strpad('INaB in component inab (A_per_F)');
    LEGEND_ALGEBRAIC(:,112) = strpad('INaB_jn in component inab (A_per_F)');
    LEGEND_ALGEBRAIC(:,117) = strpad('INaB_sl in component inab (A_per_F)');
    LEGEND_CONSTANTS(:,54) = strpad('gNaB in component inab (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,77) = strpad('INaCa in component inaca (A_per_F)');
    LEGEND_CONSTANTS(:,105) = strpad('IbarNCX in component inaca (A_per_F)');
    LEGEND_ALGEBRAIC(:,66) = strpad('Ka_jn in component inaca (dimensionless)');
    LEGEND_ALGEBRAIC(:,67) = strpad('Ka_sl in component inaca (dimensionless)');
    LEGEND_CONSTANTS(:,55) = strpad('Kdact in component inaca (mM)');
    LEGEND_CONSTANTS(:,56) = strpad('KmCai in component inaca (mM)');
    LEGEND_CONSTANTS(:,57) = strpad('KmCao in component inaca (mM)');
    LEGEND_CONSTANTS(:,58) = strpad('KmNai in component inaca (mM)');
    LEGEND_CONSTANTS(:,59) = strpad('KmNao in component inaca (mM)');
    LEGEND_CONSTANTS(:,60) = strpad('Q10NCX in component inaca (dimensionless)');
    LEGEND_ALGEBRAIC(:,68) = strpad('inaca_INaCa_jn_s1 in component inaca (mol4_per_m12)');
    LEGEND_ALGEBRAIC(:,69) = strpad('inaca_INaCa_jn_s2 in component inaca (mol4_per_m12)');
    LEGEND_ALGEBRAIC(:,70) = strpad('inaca_INaCa_jn_s3 in component inaca (mol4_per_m12)');
    LEGEND_ALGEBRAIC(:,72) = strpad('inaca_INaCa_sl_s1 in component inaca (mol4_per_m12)');
    LEGEND_ALGEBRAIC(:,73) = strpad('inaca_INaCa_sl_s2 in component inaca (mol4_per_m12)');
    LEGEND_ALGEBRAIC(:,74) = strpad('inaca_INaCa_sl_s3 in component inaca (mol4_per_m12)');
    LEGEND_CONSTANTS(:,61) = strpad('ksat in component inaca (dimensionless)');
    LEGEND_CONSTANTS(:,62) = strpad('nu in component inaca (dimensionless)');
    LEGEND_ALGEBRAIC(:,80) = strpad('INaK in component inak (A_per_F)');
    LEGEND_ALGEBRAIC(:,78) = strpad('INaK_jn in component inak (A_per_F)');
    LEGEND_ALGEBRAIC(:,79) = strpad('INaK_sl in component inak (A_per_F)');
    LEGEND_CONSTANTS(:,63) = strpad('IbarNaK in component inak (A_per_F)');
    LEGEND_CONSTANTS(:,64) = strpad('KmKo in component inak (mM)');
    LEGEND_CONSTANTS(:,65) = strpad('KmNaip in component inak (mM)');
    LEGEND_ALGEBRAIC(:,76) = strpad('fnak in component inak (dimensionless)');
    LEGEND_CONSTANTS(:,106) = strpad('sigma in component inak (dimensionless)');
    LEGEND_ALGEBRAIC(:,120) = strpad('INaL in component inal (A_per_F)');
    LEGEND_ALGEBRAIC(:,113) = strpad('INaL_jn in component inal (A_per_F)');
    LEGEND_ALGEBRAIC(:,119) = strpad('INaL_sl in component inal (A_per_F)');
    LEGEND_CONSTANTS(:,107) = strpad('gNaL in component inal (mS_per_uF)');
    LEGEND_STATES(:,35) = strpad('hl in component inal (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,10) = strpad('inal_hl_inf in component inal (dimensionless)');
    LEGEND_CONSTANTS(:,66) = strpad('inal_hl_tau in component inal (ms)');
    LEGEND_ALGEBRAIC(:,11) = strpad('inal_ml_a in component inal (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,24) = strpad('inal_ml_b in component inal (mS_per_uF)');
    LEGEND_STATES(:,36) = strpad('ml in component inal (mS_per_uF)');
    LEGEND_CONSTANTS(:,67) = strpad('Cl_i in component ion (mM)');
    LEGEND_CONSTANTS(:,68) = strpad('Cl_o in component ion (mM)');
    LEGEND_CONSTANTS(:,69) = strpad('IbarSLCaP in component ipca (A_per_F)');
    LEGEND_ALGEBRAIC(:,86) = strpad('IpCa in component ipca (A_per_F)');
    LEGEND_CONSTANTS(:,70) = strpad('KmPCa in component ipca (mM)');
    LEGEND_CONSTANTS(:,71) = strpad('Q10SLCaP in component ipca (dimensionless)');
    LEGEND_CONSTANTS(:,90) = strpad('b in component ipca (dimensionless)');
    LEGEND_ALGEBRAIC(:,81) = strpad('ipca_IpCa_jn_a in component ipca (dimensionless)');
    LEGEND_ALGEBRAIC(:,83) = strpad('ipca_IpCa_sl_a in component ipca (dimensionless)');
    LEGEND_ALGEBRAIC(:,109) = strpad('Ito in component ito (A_per_F)');
    LEGEND_CONSTANTS(:,108) = strpad('gto in component ito (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,12) = strpad('ito_x_inf in component ito (dimensionless)');
    LEGEND_ALGEBRAIC(:,25) = strpad('ito_x_tau in component ito (ms)');
    LEGEND_ALGEBRAIC(:,13) = strpad('ito_y_inf in component ito (dimensionless)');
    LEGEND_ALGEBRAIC(:,26) = strpad('ito_y_tau in component ito (ms)');
    LEGEND_STATES(:,37) = strpad('x in component ito (mS_per_uF)');
    LEGEND_STATES(:,38) = strpad('y in component ito (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,97) = strpad('ICa_tot in component membrane (A_per_F)');
    LEGEND_ALGEBRAIC(:,85) = strpad('ICl_tot in component membrane (A_per_F)');
    LEGEND_ALGEBRAIC(:,121) = strpad('IK_tot in component potassium (A_per_F)');
    LEGEND_ALGEBRAIC(:,124) = strpad('INa_tot in component membrane (A_per_F)');
    LEGEND_ALGEBRAIC(:,122) = strpad('INa_tot_jn in component sodium (A_per_F)');
    LEGEND_ALGEBRAIC(:,123) = strpad('INa_tot_sl in component sodium (A_per_F)');
    LEGEND_ALGEBRAIC(:,126) = strpad('I_stim in component stimulus (A_per_F)');
    LEGEND_ALGEBRAIC(:,125) = strpad('I_tot in component membrane (A_per_F)');
    LEGEND_CONSTANTS(:,72) = strpad('R in component phys (mJ_per_mol_per_K)');
    LEGEND_CONSTANTS(:,73) = strpad('T in component phys (kelvin)');
    LEGEND_CONSTANTS(:,109) = strpad('Kmf in component ryr (mM)');
    LEGEND_CONSTANTS(:,74) = strpad('Kmr in component ryr (mM)');
    LEGEND_CONSTANTS(:,75) = strpad('MaxSR in component ryr (dimensionless)');
    LEGEND_CONSTANTS(:,76) = strpad('MinSR in component ryr (dimensionless)');
    LEGEND_CONSTANTS(:,77) = strpad('Q10SRCaP in component ryr (dimensionless)');
    LEGEND_ALGEBRAIC(:,14) = strpad('RI in component ryr (dimensionless)');
    LEGEND_CONSTANTS(:,78) = strpad('Vmax_SRCaP in component ryr (mM_per_ms)');
    LEGEND_CONSTANTS(:,79) = strpad('ec50SR in component ryr (mM)');
    LEGEND_CONSTANTS(:,80) = strpad('hillSRCaP in component ryr (dimensionless)');
    LEGEND_STATES(:,39) = strpad('i in component ryr (mS_per_uF)');
    LEGEND_ALGEBRAIC(:,27) = strpad('kCaSR in component ryr (dimensionless)');
    LEGEND_CONSTANTS(:,81) = strpad('kiCa in component ryr (per_mM_per_ms)');
    LEGEND_ALGEBRAIC(:,31) = strpad('kiSRCa in component ryr (per_mM_per_ms)');
    LEGEND_CONSTANTS(:,82) = strpad('kim in component ryr (mS_per_uF)');
    LEGEND_CONSTANTS(:,110) = strpad('koCa in component ryr (per_mM2_per_ms)');
    LEGEND_ALGEBRAIC(:,35) = strpad('koSRCa in component ryr (per_mM2_per_ms)');
    LEGEND_CONSTANTS(:,83) = strpad('kom in component ryr (mS_per_uF)');
    LEGEND_CONSTANTS(:,84) = strpad('ks in component ryr (mS_per_uF)');
    LEGEND_STATES(:,40) = strpad('o in component ryr (mS_per_uF)');
    LEGEND_STATES(:,41) = strpad('r in component ryr (mS_per_uF)');
    LEGEND_CONSTANTS(:,111) = strpad('JNa_jnsl in component sodium (m3_per_s)');
    LEGEND_CONSTANTS(:,112) = strpad('JNa_slmyo in component sodium (m3_per_s)');
    LEGEND_CONSTANTS(:,113) = strpad('amplitude in component stimulus (A_per_F)');
    LEGEND_CONSTANTS(:,85) = strpad('duration in component stimulus (ms)');
    LEGEND_CONSTANTS(:,86) = strpad('offset in component stimulus (ms)');
    LEGEND_CONSTANTS(:,87) = strpad('period in component stimulus (ms)');
    LEGEND_RATES(:,1) = strpad('d/dt CaM in component buffca (mM)');
    LEGEND_RATES(:,14) = strpad('d/dt TnCL in component buffca (mM)');
    LEGEND_RATES(:,12) = strpad('d/dt TnCHc in component buffca (mM)');
    LEGEND_RATES(:,13) = strpad('d/dt TnCHm in component buffca (mM)');
    LEGEND_RATES(:,5) = strpad('d/dt Myoc in component buffca (mM)');
    LEGEND_RATES(:,6) = strpad('d/dt Myom in component buffca (mM)');
    LEGEND_RATES(:,11) = strpad('d/dt SRB in component buffca (mM)');
    LEGEND_RATES(:,9) = strpad('d/dt SLL_jn in component buffca (mM)');
    LEGEND_RATES(:,7) = strpad('d/dt SLH_jn in component buffca (mM)');
    LEGEND_RATES(:,10) = strpad('d/dt SLL_sl in component buffca (mM)');
    LEGEND_RATES(:,8) = strpad('d/dt SLH_sl in component buffca (mM)');
    LEGEND_RATES(:,15) = strpad('d/dt NaB_jn in component buffna (mM)');
    LEGEND_RATES(:,16) = strpad('d/dt NaB_sl in component buffna (mM)');
    LEGEND_RATES(:,2) = strpad('d/dt Ca_i in component calcium (mM)');
    LEGEND_RATES(:,3) = strpad('d/dt Ca_jn in component calcium (mM)');
    LEGEND_RATES(:,4) = strpad('d/dt Ca_sl in component calcium (mM)');
    LEGEND_RATES(:,19) = strpad('d/dt Ca_sr in component calcium (mM)');
    LEGEND_RATES(:,20) = strpad('d/dt Csqn in component calcium (mM)');
    LEGEND_RATES(:,23) = strpad('d/dt d in component ical (mS_per_uF)');
    LEGEND_RATES(:,24) = strpad('d/dt f in component ical (mS_per_uF)');
    LEGEND_RATES(:,25) = strpad('d/dt fCaB_jn in component ical (mS_per_uF)');
    LEGEND_RATES(:,26) = strpad('d/dt fCaB_sl in component ical (mS_per_uF)');
    LEGEND_RATES(:,27) = strpad('d/dt x in component ikr (mS_per_uF)');
    LEGEND_RATES(:,29) = strpad('d/dt x in component iks (mS_per_uF)');
    LEGEND_RATES(:,30) = strpad('d/dt r in component ikur (mS_per_uF)');
    LEGEND_RATES(:,31) = strpad('d/dt s in component ikur (mS_per_uF)');
    LEGEND_RATES(:,32) = strpad('d/dt h in component ina (mS_per_uF)');
    LEGEND_RATES(:,33) = strpad('d/dt j in component ina (mS_per_uF)');
    LEGEND_RATES(:,34) = strpad('d/dt m in component ina (mS_per_uF)');
    LEGEND_RATES(:,35) = strpad('d/dt hl in component inal (mS_per_uF)');
    LEGEND_RATES(:,36) = strpad('d/dt ml in component inal (mS_per_uF)');
    LEGEND_RATES(:,37) = strpad('d/dt x in component ito (mS_per_uF)');
    LEGEND_RATES(:,38) = strpad('d/dt y in component ito (mS_per_uF)');
    LEGEND_RATES(:,21) = strpad('d/dt V in component membrane (mV)');
    LEGEND_RATES(:,22) = strpad('d/dt K_i in component potassium (mM)');
    LEGEND_RATES(:,39) = strpad('d/dt i in component ryr (mS_per_uF)');
    LEGEND_RATES(:,40) = strpad('d/dt o in component ryr (mS_per_uF)');
    LEGEND_RATES(:,41) = strpad('d/dt r in component ryr (mS_per_uF)');
    LEGEND_RATES(:,28) = strpad('d/dt Na_i in component sodium (mM)');
    LEGEND_RATES(:,17) = strpad('d/dt Na_jn in component sodium (mM)');
    LEGEND_RATES(:,18) = strpad('d/dt Na_sl in component sodium (mM)');
    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) = 0.024;
    CONSTANTS(:,2) = 0.14;
    CONSTANTS(:,3) = 0.07;
    CONSTANTS(:,4) = 0.14;
    STATES(:,1) = 6.73710089221736985e-04;
    STATES(:,2) = 2.02054825252250719e-04;
    STATES(:,3) = 3.07333085373605341e-04;
    STATES(:,4) = 2.23134214147860355e-04;
    CONSTANTS(:,5) = 1.0;
    STATES(:,5) = 3.79608844523014136e-03;
    STATES(:,6) = 1.35691489581683838e-01;
    STATES(:,7) = 1.00754933940878755e-01;
    STATES(:,8) = 1.86055242143676608e-01;
    STATES(:,9) = 1.28122161143722785e-02;
    STATES(:,10) = 2.05131045874987414e-02;
    STATES(:,11) = 4.31383413366585346e-03;
    STATES(:,12) = 1.27405935029765827e-01;
    STATES(:,13) = 5.91106798633337510e-03;
    STATES(:,14) = 1.77305527980331397e-02;
    CONSTANTS(:,6) = 0.238;
    CONSTANTS(:,7) = 0.00046;
    CONSTANTS(:,8) = 5.7e-05;
    CONSTANTS(:,9) = 0.03;
    CONSTANTS(:,10) = 1.3;
    CONSTANTS(:,11) = 0.06;
    CONSTANTS(:,12) = 3.2e-05;
    CONSTANTS(:,13) = 0.00333;
    CONSTANTS(:,14) = 0.0196;
    CONSTANTS(:,15) = 34.0;
    CONSTANTS(:,16) = 13.8;
    CONSTANTS(:,17) = 0.0157;
    CONSTANTS(:,18) = 100.0;
    CONSTANTS(:,19) = 100.0;
    CONSTANTS(:,20) = 100.0;
    CONSTANTS(:,21) = 2.37;
    CONSTANTS(:,22) = 0.003;
    CONSTANTS(:,23) = 32.7;
    CONSTANTS(:,24) = 7.561;
    CONSTANTS(:,25) = 1.65;
    STATES(:,15) = 3.62527560419928152e+00;
    STATES(:,16) = 7.91087217252140484e-01;
    STATES(:,17) = 9.20484970072274145e+00;
    STATES(:,18) = 9.20559631626608166e+00;
    CONSTANTS(:,26) = 0.001;
    CONSTANTS(:,27) = 0.0001;
    CONSTANTS(:,28) = 1.381e-10;
    STATES(:,19) = 4.85784491491899340e-01;
    STATES(:,20) = 1.11204132576763315e+00;
    CONSTANTS(:,29) = 96485.0;
    CONSTANTS(:,30) = 65.0;
    CONSTANTS(:,31) = 100.0;
    CONSTANTS(:,32) = 0.0;
    CONSTANTS(:,33) = 100.0;
    CONSTANTS(:,34) = 10.25;
    CONSTANTS(:,35) = 3.14159265358979312e+00;
    CONSTANTS(:,36) = 0.11;
    STATES(:,21) = -7.49738877185916408e+01;
    CONSTANTS(:,37) = 6.06430000000000033e-04;
    CONSTANTS(:,38) = 1.8;
    CONSTANTS(:,39) = 0.9;
    STATES(:,22) = 120.0;
    CONSTANTS(:,40) = 5.4;
    CONSTANTS(:,41) = 140.0;
    CONSTANTS(:,42) = 1.8;
    STATES(:,23) = 1.67750282273476489e-05;
    STATES(:,24) = 9.98760022198537767e-01;
    STATES(:,25) = 4.24066845471226944e-02;
    STATES(:,26) = 3.13717989771996886e-02;
    CONSTANTS(:,43) = 0.0;
    CONSTANTS(:,44) = 0.0;
    CONSTANTS(:,45) = 0.009;
    CONSTANTS(:,46) = 0.1;
    CONSTANTS(:,47) = 0.0548;
    CONSTANTS(:,48) = 2.1;
    CONSTANTS(:,49) = 0.0;
    CONSTANTS(:,50) = 0.0;
    CONSTANTS(:,51) = 0.002;
    STATES(:,27) = 4.42405976310802831e-03;
    STATES(:,28) = 9.20582356442186978e+00;
    CONSTANTS(:,52) = 0.01833;
    STATES(:,29) = 6.73017121702351208e-03;
    STATES(:,30) = 3.28872423495696863e-04;
    STATES(:,31) = 9.62950458643225149e-01;
    CONSTANTS(:,53) = 23.0;
    STATES(:,32) = 8.70955463748630310e-01;
    STATES(:,33) = 8.92501441492474679e-01;
    STATES(:,34) = 7.97426794491632243e-03;
    CONSTANTS(:,54) = 0.000597;
    CONSTANTS(:,55) = 0.000384;
    CONSTANTS(:,56) = 0.00359;
    CONSTANTS(:,57) = 1.3;
    CONSTANTS(:,58) = 12.29;
    CONSTANTS(:,59) = 87.5;
    CONSTANTS(:,60) = 1.57;
    CONSTANTS(:,61) = 0.27;
    CONSTANTS(:,62) = 0.35;
    CONSTANTS(:,63) = 1.26;
    CONSTANTS(:,64) = 1.5;
    CONSTANTS(:,65) = 11.0;
    STATES(:,35) = 5.18117114324953484e-02;
    CONSTANTS(:,66) = 600.0;
    STATES(:,36) = 7.97426794490581868e-03;
    CONSTANTS(:,67) = 15.0;
    CONSTANTS(:,68) = 150.0;
    CONSTANTS(:,69) = 0.0471;
    CONSTANTS(:,70) = 0.0005;
    CONSTANTS(:,71) = 2.35;
    STATES(:,37) = 1.19930453889315167e-03;
    STATES(:,38) = 9.52396284498915668e-01;
    CONSTANTS(:,72) = 8314.0;
    CONSTANTS(:,73) = 310.0;
    CONSTANTS(:,74) = 1.7;
    CONSTANTS(:,75) = 15.0;
    CONSTANTS(:,76) = 1.0;
    CONSTANTS(:,77) = 2.6;
    CONSTANTS(:,78) = 0.0053114;
    CONSTANTS(:,79) = 0.45;
    CONSTANTS(:,80) = 1.787;
    STATES(:,39) = 4.16351235655844508e-07;
    CONSTANTS(:,81) = 0.5;
    CONSTANTS(:,82) = 0.005;
    CONSTANTS(:,83) = 0.06;
    CONSTANTS(:,84) = 25.0;
    STATES(:,40) = 1.73377710603823675e-06;
    STATES(:,41) = 8.06359263514647706e-01;
    CONSTANTS(:,85) = 0.5;
    CONSTANTS(:,86) = 50.0;
    CONSTANTS(:,87) = 1000.0;
    CONSTANTS(:,88) =  19.0000.*0.000900000;
    CONSTANTS(:,89) =  ( ( CONSTANTS(:,35).*power(CONSTANTS(:,34), 2.00000)).*CONSTANTS(:,33)).*1.00000e-15;
    CONSTANTS(:,90) = power( CONSTANTS(:,70).*1.00000, 1.60000);
    CONSTANTS(:,91) = 1.00000./1.21340e+12;
    CONSTANTS(:,92) = 1.00000./2.68510e+11;
    CONSTANTS(:,93) = (CONSTANTS(:,29)./CONSTANTS(:,72))./CONSTANTS(:,73);
    CONSTANTS(:,94) = 1.00000 - CONSTANTS(:,39);
    CONSTANTS(:,95) =  (1.00000 -  0.500000.*CONSTANTS(:,32)).*0.000270000;
    CONSTANTS(:,96) =  (1.00000 -  0.500000.*CONSTANTS(:,32)).*1.35000e-07;
    CONSTANTS(:,97) =  (1.00000 -  0.500000.*CONSTANTS(:,32)).*7.50000e-09;
    CONSTANTS(:,98) =  ( (1.00000+CONSTANTS(:,32)).*power((CONSTANTS(:,40)./5.40000), 1.0 ./ 2)).*0.0525000;
    CONSTANTS(:,99) = CONSTANTS(:,50)+( (1.00000 - CONSTANTS(:,50)).*CONSTANTS(:,49))./(CONSTANTS(:,49)+0.125000);
    CONSTANTS(:,100) = 0.100000 -  0.0500000.*CONSTANTS(:,32);
    CONSTANTS(:,101) =  0.0350000.*power((CONSTANTS(:,40)./5.40000), 1.0 ./ 2);
    CONSTANTS(:,102) =  (1.00000+CONSTANTS(:,32)).*0.00350000;
    CONSTANTS(:,103) =  (1.00000+CONSTANTS(:,32)).*0.00350000;
    CONSTANTS(:,104) =  (1.00000 -  0.500000.*CONSTANTS(:,32)).*0.0450000;
    CONSTANTS(:,105) =  (1.00000+ 0.400000.*CONSTANTS(:,32)).*3.15000;
    CONSTANTS(:,106) = (exp(CONSTANTS(:,41)./67.3000) - 1.00000)./7.00000;
    CONSTANTS(:,107) =  0.00250000.*CONSTANTS(:,32);
    CONSTANTS(:,108) =  (1.00000 -  0.700000.*CONSTANTS(:,32)).*0.165000;
    CONSTANTS(:,109) =  2.50000.*0.000246000;
    CONSTANTS(:,110) =  (10.0000+ 20.0000.*CONSTANTS(:,32)).*1.00000;
    CONSTANTS(:,111) = 1.00000./( (1.63820e+12./3.00000).*100.000);
    CONSTANTS(:,112) = 1.00000./( (1.83080e+10./3.00000).*100.000);
    CONSTANTS(:,113) =  - 90.0000;
    CONSTANTS(:,125) = 0.00000;
    CONSTANTS(:,114) =  ( 0.0539000.*0.0100000).*CONSTANTS(:,89);
    CONSTANTS(:,115) = 1.00000 - CONSTANTS(:,36);
    CONSTANTS(:,116) =  (1.00000./CONSTANTS(:,93)).*log(CONSTANTS(:,67)./CONSTANTS(:,68));
    CONSTANTS(:,117) =  0.650000.*CONSTANTS(:,89);
    CONSTANTS(:,118) = (CONSTANTS(:,73) - 310.000)./10.0000;
    CONSTANTS(:,119) =  (( 0.00165000.*CONSTANTS(:,117))./CONSTANTS(:,114)).*0.100000;
    CONSTANTS(:,120) =  (( 0.00460000.*CONSTANTS(:,117))./CONSTANTS(:,114)).*0.100000;
    CONSTANTS(:,121) =  0.0200000.*CONSTANTS(:,89);
    CONSTANTS(:,122) =  0.0350000.*CONSTANTS(:,89);
    CONSTANTS(:,123) =  (CONSTANTS(:,117)./CONSTANTS(:,121)).*0.0134000;
    CONSTANTS(:,124) =  (CONSTANTS(:,117)./CONSTANTS(:,121)).*0.0374000;
    CONSTANTS(:,125) =  (CONSTANTS(:,117)./CONSTANTS(:,122)).*0.140000;
    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(:,22) = CONSTANTS(:,125);
    RATES(:,1) =  ( CONSTANTS(:,15).*STATES(:,2)).*(CONSTANTS(:,1) - STATES(:,1)) -  CONSTANTS(:,6).*STATES(:,1);
    RATES(:,25) =  ( ( 1.70000.*STATES(:,3)).*(1.00000 - STATES(:,25)) -  0.0119000.*STATES(:,25)).*1.00000;
    RATES(:,26) =  ( ( 1.70000.*STATES(:,4)).*(1.00000 - STATES(:,26)) -  0.0119000.*STATES(:,26)).*1.00000;
    RATES(:,28) =  (CONSTANTS(:,112)./CONSTANTS(:,117)).*(STATES(:,18) - STATES(:,28));
    RATES(:,5) =  ( CONSTANTS(:,16).*STATES(:,2)).*((CONSTANTS(:,4) - STATES(:,5)) - STATES(:,6)) -  CONSTANTS(:,7).*STATES(:,5);
    ALGEBRAIC(:,10) = 1.00000./(1.00000+exp((STATES(:,21)+91.0000)./6.10000));
    RATES(:,35) = (ALGEBRAIC(:,10) - STATES(:,35))./CONSTANTS(:,66);
    RATES(:,6) =  ( CONSTANTS(:,17).*CONSTANTS(:,5)).*((CONSTANTS(:,4) - STATES(:,5)) - STATES(:,6)) -  CONSTANTS(:,8).*STATES(:,6);
    ALGEBRAIC(:,1) = 1.00000./(1.00000+exp( - (STATES(:,21)+9.00000)./6.00000));
    ALGEBRAIC(:,15) = ( ALGEBRAIC(:,1).*(1.00000 - exp( - (STATES(:,21)+9.00000)./6.00000)))./( 0.0350000.*(STATES(:,21)+9.00000));
    RATES(:,23) = (ALGEBRAIC(:,1) - STATES(:,23))./ALGEBRAIC(:,15);
    ALGEBRAIC(:,2) = 1.00000./(1.00000+exp((STATES(:,21)+30.0000)./7.00000))+0.200000./(1.00000+exp((50.0000 - STATES(:,21))./20.0000));
    ALGEBRAIC(:,16) = 1.00000./( 0.0197000.*exp( - power( 0.0337000.*(STATES(:,21)+25.0000), 2.00000))+0.0200000);
    RATES(:,24) = (ALGEBRAIC(:,2) - STATES(:,24))./ALGEBRAIC(:,16);
    ALGEBRAIC(:,3) = 1.00000./(1.00000+exp( - (STATES(:,21)+10.0000)./5.00000));
    ALGEBRAIC(:,17) = ( ((550.000)./(1.00000+exp(( - 22.0000 - STATES(:,21))./9.00000))).*6.00000)./(1.00000+exp((STATES(:,21)+11.0000)./9.00000))+230.000./(1.00000+exp((STATES(:,21)+40.0000)./20.0000));
    RATES(:,27) = (ALGEBRAIC(:,3) - STATES(:,27))./ALGEBRAIC(:,17);
    ALGEBRAIC(:,4) = 1.00000./(1.00000+exp( - (STATES(:,21)+3.80000)./14.2500));
    ALGEBRAIC(:,18) = 990.100./(1.00000+exp( - (STATES(:,21)+2.43600)./14.1200));
    RATES(:,29) = (ALGEBRAIC(:,4) - STATES(:,29))./ALGEBRAIC(:,18);
    ALGEBRAIC(:,5) = 1.00000./(1.00000+exp((STATES(:,21)+6.00000)./ - 8.60000));
    ALGEBRAIC(:,19) = 9.00000./(1.00000+exp((STATES(:,21)+5.00000)./12.0000))+0.500000;
    RATES(:,30) = (ALGEBRAIC(:,5) - STATES(:,30))./ALGEBRAIC(:,19);
    ALGEBRAIC(:,6) = 1.00000./(1.00000+exp((STATES(:,21)+7.50000)./10.0000));
    ALGEBRAIC(:,20) = 590.000./(1.00000+exp((STATES(:,21)+60.0000)./10.0000))+3050.00;
    RATES(:,31) = (ALGEBRAIC(:,6) - STATES(:,31))./ALGEBRAIC(:,20);
    ALGEBRAIC(:,11) = ( 0.320000.*(STATES(:,21)+47.1300))./(1.00000 - exp(  - 0.100000.*(STATES(:,21)+47.1300)));
    ALGEBRAIC(:,24) =  0.0800000.*exp( - STATES(:,21)./11.0000);
    RATES(:,36) =  ALGEBRAIC(:,11).*(1.00000 - STATES(:,36)) -  ALGEBRAIC(:,24).*STATES(:,36);
    ALGEBRAIC(:,12) = 1.00000./(1.00000+exp( - (STATES(:,21)+1.00000)./11.0000));
    ALGEBRAIC(:,25) =  3.50000.*exp( - power(STATES(:,21)./30.0000, 2.00000))+1.50000;
    RATES(:,37) = (ALGEBRAIC(:,12) - STATES(:,37))./ALGEBRAIC(:,25);
    ALGEBRAIC(:,13) = 1.00000./(1.00000+exp((STATES(:,21)+40.5000)./11.5000));
    ALGEBRAIC(:,26) =  25.6350.*exp( - power((STATES(:,21)+52.4500)./15.8827, 2.00000))+24.1400;
    RATES(:,38) = (ALGEBRAIC(:,13) - STATES(:,38))./ALGEBRAIC(:,26);
    RATES(:,7) =  ( CONSTANTS(:,18).*STATES(:,3)).*(CONSTANTS(:,119) - STATES(:,7)) -  CONSTANTS(:,9).*STATES(:,7);
    RATES(:,8) =  ( CONSTANTS(:,18).*STATES(:,4)).*(CONSTANTS(:,123) - STATES(:,8)) -  CONSTANTS(:,9).*STATES(:,8);
    ALGEBRAIC(:,7) = piecewise({STATES(:,21)< - 40.0000,  0.135000.*exp((STATES(:,21)+80.0000)./ - 6.80000) }, 0.00000);
    ALGEBRAIC(:,21) = piecewise({STATES(:,21)< - 40.0000,  3.56000.*exp( 0.0790000.*STATES(:,21))+ 310000..*exp( 0.350000.*STATES(:,21)) }, 1.00000./( 0.130000.*(1.00000+exp((STATES(:,21)+10.6600)./ - 11.1000))));
    ALGEBRAIC(:,28) = ALGEBRAIC(:,7)./(ALGEBRAIC(:,7)+ALGEBRAIC(:,21));
    ALGEBRAIC(:,32) = 1.00000./(ALGEBRAIC(:,7)+ALGEBRAIC(:,21));
    RATES(:,32) = (ALGEBRAIC(:,28) - STATES(:,32))./ALGEBRAIC(:,32);
    ALGEBRAIC(:,8) = piecewise({STATES(:,21)< - 40.0000, ( (  - 127140..*exp( 0.244400.*STATES(:,21)) -  3.47400e-05.*exp(  - 0.0439100.*STATES(:,21))).*(STATES(:,21)+37.7800))./(1.00000+exp( 0.311000.*(STATES(:,21)+79.2300))) }, 0.00000);
    ALGEBRAIC(:,22) = piecewise({STATES(:,21)< - 40.0000, ( 0.121200.*exp(  - 0.0105200.*STATES(:,21)))./(1.00000+exp(  - 0.137800.*(STATES(:,21)+40.1400))) }, ( 0.300000.*exp(  - 2.53500e-07.*STATES(:,21)))./(1.00000+exp(  - 0.100000.*(STATES(:,21)+32.0000))));
    ALGEBRAIC(:,29) = ALGEBRAIC(:,8)./(ALGEBRAIC(:,8)+ALGEBRAIC(:,22));
    ALGEBRAIC(:,33) = 1.00000./(ALGEBRAIC(:,8)+ALGEBRAIC(:,22));
    RATES(:,33) = (ALGEBRAIC(:,29) - STATES(:,33))./ALGEBRAIC(:,33);
    ALGEBRAIC(:,9) = piecewise({STATES(:,21)== - 47.1300, 3.20000 }, ( 0.320000.*(STATES(:,21)+47.1300))./(1.00000 - exp(  - 0.100000.*(STATES(:,21)+47.1300))));
    ALGEBRAIC(:,23) =  0.0800000.*exp( - STATES(:,21)./11.0000);
    ALGEBRAIC(:,30) = ALGEBRAIC(:,9)./(ALGEBRAIC(:,9)+ALGEBRAIC(:,23));
    ALGEBRAIC(:,34) = 1.00000./(ALGEBRAIC(:,9)+ALGEBRAIC(:,23));
    RATES(:,34) = (ALGEBRAIC(:,30) - STATES(:,34))./ALGEBRAIC(:,34);
    ALGEBRAIC(:,14) = ((1.00000 - STATES(:,41)) - STATES(:,40)) - STATES(:,39);
    ALGEBRAIC(:,27) = CONSTANTS(:,75) - (CONSTANTS(:,75) - CONSTANTS(:,76))./(1.00000+power(CONSTANTS(:,79)./STATES(:,19), 2.50000));
    ALGEBRAIC(:,31) =  CONSTANTS(:,81).*ALGEBRAIC(:,27);
    ALGEBRAIC(:,35) = CONSTANTS(:,110)./ALGEBRAIC(:,27);
    RATES(:,39) = ( ( ALGEBRAIC(:,31).*STATES(:,3)).*STATES(:,40) -  CONSTANTS(:,82).*STATES(:,39)) - ( CONSTANTS(:,83).*STATES(:,39) -  ( ALGEBRAIC(:,35).*power(STATES(:,3), 2.00000)).*ALGEBRAIC(:,14));
    RATES(:,40) = ( ( ALGEBRAIC(:,35).*power(STATES(:,3), 2.00000)).*STATES(:,41) -  CONSTANTS(:,83).*STATES(:,40)) - ( ( ALGEBRAIC(:,31).*STATES(:,3)).*STATES(:,40) -  CONSTANTS(:,82).*STATES(:,39));
    RATES(:,41) = ( CONSTANTS(:,82).*ALGEBRAIC(:,14) -  ( ALGEBRAIC(:,31).*STATES(:,3)).*STATES(:,41)) - ( ( ALGEBRAIC(:,35).*power(STATES(:,3), 2.00000)).*STATES(:,41) -  CONSTANTS(:,83).*STATES(:,40));
    RATES(:,9) =  ( CONSTANTS(:,19).*STATES(:,3)).*(CONSTANTS(:,120) - STATES(:,9)) -  CONSTANTS(:,10).*STATES(:,9);
    RATES(:,10) =  ( CONSTANTS(:,19).*STATES(:,4)).*(CONSTANTS(:,124) - STATES(:,10)) -  CONSTANTS(:,10).*STATES(:,10);
    RATES(:,11) =  ( CONSTANTS(:,20).*STATES(:,2)).*(CONSTANTS(:,88) - STATES(:,11)) -  CONSTANTS(:,11).*STATES(:,11);
    RATES(:,12) =  ( CONSTANTS(:,21).*STATES(:,2)).*((CONSTANTS(:,2) - STATES(:,12)) - STATES(:,13)) -  CONSTANTS(:,12).*STATES(:,12);
    RATES(:,13) =  ( CONSTANTS(:,22).*CONSTANTS(:,5)).*((CONSTANTS(:,2) - STATES(:,12)) - STATES(:,13)) -  CONSTANTS(:,13).*STATES(:,13);
    RATES(:,14) =  ( CONSTANTS(:,23).*STATES(:,2)).*(CONSTANTS(:,3) - STATES(:,14)) -  CONSTANTS(:,14).*STATES(:,14);
    RATES(:,15) =  ( CONSTANTS(:,27).*STATES(:,17)).*(CONSTANTS(:,24) - STATES(:,15)) -  CONSTANTS(:,26).*STATES(:,15);
    RATES(:,16) =  ( CONSTANTS(:,27).*STATES(:,18)).*(CONSTANTS(:,25) - STATES(:,16)) -  CONSTANTS(:,26).*STATES(:,16);
    RATES(:,20) =  ( CONSTANTS(:,31).*STATES(:,19)).*(CONSTANTS(:,125) - STATES(:,20)) -  CONSTANTS(:,30).*STATES(:,20);
    ALGEBRAIC(:,36) = RATES(:,9)+RATES(:,7);
    ALGEBRAIC(:,87) =  ((1.00000./CONSTANTS(:,93))./2.00000).*log(CONSTANTS(:,38)./STATES(:,3));
    ALGEBRAIC(:,88) =  ( CONSTANTS(:,36).*CONSTANTS(:,37)).*(STATES(:,21) - ALGEBRAIC(:,87));
    ALGEBRAIC(:,39) = ( ( ( ( ( CONSTANTS(:,95).*4.00000).*STATES(:,21)).*CONSTANTS(:,29)).*CONSTANTS(:,93)).*( ( 0.341000.*STATES(:,3)).*exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,93)) -  0.341000.*CONSTANTS(:,38)))./(exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,93)) - 1.00000);
    ALGEBRAIC(:,40) =  ( ( ( ( ( CONSTANTS(:,39).*ALGEBRAIC(:,39)).*STATES(:,23)).*STATES(:,24)).*((1.00000 - STATES(:,25))+CONSTANTS(:,44))).*power(CONSTANTS(:,42), CONSTANTS(:,118))).*0.450000;
    ALGEBRAIC(:,66) = 1.00000./(1.00000+power(CONSTANTS(:,55)./STATES(:,3), 2.00000));
    ALGEBRAIC(:,68) =  ( exp( ( CONSTANTS(:,62).*STATES(:,21)).*CONSTANTS(:,93)).*power(STATES(:,17), 3.00000)).*CONSTANTS(:,38);
    ALGEBRAIC(:,69) =  ( exp( ( (CONSTANTS(:,62) - 1.00000).*STATES(:,21)).*CONSTANTS(:,93)).*power(CONSTANTS(:,41), 3.00000)).*STATES(:,3);
    ALGEBRAIC(:,70) = ((( ( CONSTANTS(:,56).*power(CONSTANTS(:,41), 3.00000)).*(1.00000+power(STATES(:,17)./CONSTANTS(:,58), 3.00000))+ ( power(CONSTANTS(:,59), 3.00000).*STATES(:,3)).*(1.00000+STATES(:,3)./CONSTANTS(:,56)))+ CONSTANTS(:,57).*power(STATES(:,17), 3.00000))+ power(STATES(:,17), 3.00000).*CONSTANTS(:,38))+ power(CONSTANTS(:,41), 3.00000).*STATES(:,3);
    ALGEBRAIC(:,71) = (( ( ( ( CONSTANTS(:,36).*CONSTANTS(:,105)).*power(CONSTANTS(:,60), CONSTANTS(:,118))).*ALGEBRAIC(:,66)).*(ALGEBRAIC(:,68) - ALGEBRAIC(:,69)))./ALGEBRAIC(:,70))./(1.00000+ CONSTANTS(:,61).*exp( ( (CONSTANTS(:,62) - 1.00000).*STATES(:,21)).*CONSTANTS(:,93)));
    ALGEBRAIC(:,81) = power( STATES(:,3).*1.00000, 1.60000);
    ALGEBRAIC(:,82) = ( ( ( CONSTANTS(:,36).*power(CONSTANTS(:,71), CONSTANTS(:,118))).*CONSTANTS(:,69)).*ALGEBRAIC(:,81))./(ALGEBRAIC(:,81)+CONSTANTS(:,90));
    ALGEBRAIC(:,89) = ((ALGEBRAIC(:,40)+ALGEBRAIC(:,88))+ALGEBRAIC(:,82)) -  2.00000.*ALGEBRAIC(:,71);
    ALGEBRAIC(:,90) =  ( CONSTANTS(:,84).*STATES(:,40)).*(STATES(:,19) - STATES(:,3));
    ALGEBRAIC(:,92) =  ( (1.00000+ 0.250000.*CONSTANTS(:,32)).*(STATES(:,19) - STATES(:,3))).*5.34800e-06;
    RATES(:,3) = ((((  - ALGEBRAIC(:,89).*CONSTANTS(:,28))./( ( CONSTANTS(:,114).*2.00000).*CONSTANTS(:,29))+ (CONSTANTS(:,91)./CONSTANTS(:,114)).*(STATES(:,4) - STATES(:,3))) - ALGEBRAIC(:,36))+( ALGEBRAIC(:,90).*CONSTANTS(:,122))./CONSTANTS(:,114))+( ALGEBRAIC(:,92).*CONSTANTS(:,117))./CONSTANTS(:,114);
    ALGEBRAIC(:,38) = (((((RATES(:,14)+RATES(:,12))+RATES(:,13))+RATES(:,1))+RATES(:,5))+RATES(:,6))+RATES(:,11);
    ALGEBRAIC(:,94) = ( ( power(CONSTANTS(:,77), CONSTANTS(:,118)).*CONSTANTS(:,78)).*(power(STATES(:,2)./CONSTANTS(:,109), CONSTANTS(:,80)) - power(STATES(:,19)./CONSTANTS(:,74), CONSTANTS(:,80))))./((1.00000+power(STATES(:,2)./CONSTANTS(:,109), CONSTANTS(:,80)))+power(STATES(:,19)./CONSTANTS(:,74), CONSTANTS(:,80)));
    RATES(:,2) = ((  - ALGEBRAIC(:,94).*CONSTANTS(:,122))./CONSTANTS(:,117) - ALGEBRAIC(:,38))+ (CONSTANTS(:,92)./CONSTANTS(:,117)).*(STATES(:,4) - STATES(:,2));
    ALGEBRAIC(:,37) = RATES(:,10)+RATES(:,8);
    ALGEBRAIC(:,91) =  ((1.00000./CONSTANTS(:,93))./2.00000).*log(CONSTANTS(:,38)./STATES(:,4));
    ALGEBRAIC(:,93) =  ( CONSTANTS(:,115).*CONSTANTS(:,37)).*(STATES(:,21) - ALGEBRAIC(:,91));
    ALGEBRAIC(:,41) = ( ( ( ( ( CONSTANTS(:,95).*4.00000).*STATES(:,21)).*CONSTANTS(:,29)).*CONSTANTS(:,93)).*( ( 0.341000.*STATES(:,4)).*exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,93)) -  0.341000.*CONSTANTS(:,38)))./(exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,93)) - 1.00000);
    ALGEBRAIC(:,42) =  ( ( ( ( ( CONSTANTS(:,94).*ALGEBRAIC(:,41)).*STATES(:,23)).*STATES(:,24)).*((1.00000 - STATES(:,26))+CONSTANTS(:,43))).*power(CONSTANTS(:,42), CONSTANTS(:,118))).*0.450000;
    ALGEBRAIC(:,67) = 1.00000./(1.00000+power(CONSTANTS(:,55)./STATES(:,4), 2.00000));
    ALGEBRAIC(:,72) =  ( exp( ( CONSTANTS(:,62).*STATES(:,21)).*CONSTANTS(:,93)).*power(STATES(:,18), 3.00000)).*CONSTANTS(:,38);
    ALGEBRAIC(:,73) =  ( exp( ( (CONSTANTS(:,62) - 1.00000).*STATES(:,21)).*CONSTANTS(:,93)).*power(CONSTANTS(:,41), 3.00000)).*STATES(:,4);
    ALGEBRAIC(:,74) = ((( ( CONSTANTS(:,56).*power(CONSTANTS(:,41), 3.00000)).*(1.00000+power(STATES(:,18)./CONSTANTS(:,58), 3.00000))+ ( power(CONSTANTS(:,59), 3.00000).*STATES(:,4)).*(1.00000+STATES(:,4)./CONSTANTS(:,56)))+ CONSTANTS(:,57).*power(STATES(:,18), 3.00000))+ power(STATES(:,18), 3.00000).*CONSTANTS(:,38))+ power(CONSTANTS(:,41), 3.00000).*STATES(:,4);
    ALGEBRAIC(:,75) = (( ( ( ( CONSTANTS(:,115).*CONSTANTS(:,105)).*power(CONSTANTS(:,60), CONSTANTS(:,118))).*ALGEBRAIC(:,67)).*(ALGEBRAIC(:,72) - ALGEBRAIC(:,73)))./ALGEBRAIC(:,74))./(1.00000+ CONSTANTS(:,61).*exp( ( (CONSTANTS(:,62) - 1.00000).*STATES(:,21)).*CONSTANTS(:,93)));
    ALGEBRAIC(:,83) = power( STATES(:,4).*1.00000, 1.60000);
    ALGEBRAIC(:,84) = ( ( ( CONSTANTS(:,115).*power(CONSTANTS(:,71), CONSTANTS(:,118))).*CONSTANTS(:,69)).*ALGEBRAIC(:,83))./(ALGEBRAIC(:,83)+CONSTANTS(:,90));
    ALGEBRAIC(:,96) = ((ALGEBRAIC(:,42)+ALGEBRAIC(:,93))+ALGEBRAIC(:,84)) -  2.00000.*ALGEBRAIC(:,75);
    RATES(:,4) = (((  - ALGEBRAIC(:,96).*CONSTANTS(:,28))./( ( CONSTANTS(:,121).*2.00000).*CONSTANTS(:,29))+ (CONSTANTS(:,91)./CONSTANTS(:,121)).*(STATES(:,3) - STATES(:,4)))+ (CONSTANTS(:,92)./CONSTANTS(:,121)).*(STATES(:,2) - STATES(:,4))) - ALGEBRAIC(:,37);
    RATES(:,19) = (ALGEBRAIC(:,94) - (( ALGEBRAIC(:,92).*CONSTANTS(:,117))./CONSTANTS(:,122)+ALGEBRAIC(:,90))) - RATES(:,20);
    ALGEBRAIC(:,46) = ( ( ( ( CONSTANTS(:,97).*STATES(:,21)).*CONSTANTS(:,29)).*CONSTANTS(:,93)).*( ( 0.750000.*STATES(:,17)).*exp( STATES(:,21).*CONSTANTS(:,93)) -  0.750000.*CONSTANTS(:,41)))./(exp( STATES(:,21).*CONSTANTS(:,93)) - 1.00000);
    ALGEBRAIC(:,47) =  ( ( ( ( ( CONSTANTS(:,39).*ALGEBRAIC(:,46)).*STATES(:,23)).*STATES(:,24)).*((1.00000 - STATES(:,25))+CONSTANTS(:,44))).*power(CONSTANTS(:,42), CONSTANTS(:,118))).*0.450000;
    ALGEBRAIC(:,110) =  (1.00000./CONSTANTS(:,93)).*log(CONSTANTS(:,41)./STATES(:,17));
    ALGEBRAIC(:,111) =  ( ( ( ( CONSTANTS(:,36).*CONSTANTS(:,53)).*power(STATES(:,34), 3.00000)).*STATES(:,32)).*STATES(:,33)).*(STATES(:,21) - ALGEBRAIC(:,110));
    ALGEBRAIC(:,112) =  ( CONSTANTS(:,36).*CONSTANTS(:,54)).*(STATES(:,21) - ALGEBRAIC(:,110));
    ALGEBRAIC(:,76) = 1.00000./((1.00000+ 0.124500.*exp( (  - 0.100000.*STATES(:,21)).*CONSTANTS(:,93)))+ ( 0.0365000.*CONSTANTS(:,106)).*exp(  - STATES(:,21).*CONSTANTS(:,93)));
    ALGEBRAIC(:,78) = (( ( ( CONSTANTS(:,36).*CONSTANTS(:,63)).*ALGEBRAIC(:,76)).*CONSTANTS(:,40))./(1.00000+power(CONSTANTS(:,65)./STATES(:,17), 4.00000)))./(CONSTANTS(:,40)+CONSTANTS(:,64));
    ALGEBRAIC(:,113) =  ( ( ( CONSTANTS(:,36).*CONSTANTS(:,107)).*power(STATES(:,36), 3.00000)).*STATES(:,35)).*(STATES(:,21) - ALGEBRAIC(:,110));
    ALGEBRAIC(:,122) = ((((ALGEBRAIC(:,111)+ALGEBRAIC(:,112))+ 3.00000.*ALGEBRAIC(:,71))+ 3.00000.*ALGEBRAIC(:,78))+ALGEBRAIC(:,47))+ALGEBRAIC(:,113);
    RATES(:,17) = ((  - ALGEBRAIC(:,122).*CONSTANTS(:,28))./( CONSTANTS(:,114).*CONSTANTS(:,29))+ (CONSTANTS(:,111)./CONSTANTS(:,114)).*(STATES(:,18) - STATES(:,17))) - RATES(:,15);
    ALGEBRAIC(:,48) = ( ( ( ( CONSTANTS(:,97).*STATES(:,21)).*CONSTANTS(:,29)).*CONSTANTS(:,93)).*( ( 0.750000.*STATES(:,18)).*exp( STATES(:,21).*CONSTANTS(:,93)) -  0.750000.*CONSTANTS(:,41)))./(exp( STATES(:,21).*CONSTANTS(:,93)) - 1.00000);
    ALGEBRAIC(:,49) =  ( ( ( ( ( CONSTANTS(:,94).*ALGEBRAIC(:,48)).*STATES(:,23)).*STATES(:,24)).*((1.00000 - STATES(:,26))+CONSTANTS(:,43))).*power(CONSTANTS(:,42), CONSTANTS(:,118))).*0.450000;
    ALGEBRAIC(:,114) =  (1.00000./CONSTANTS(:,93)).*log(CONSTANTS(:,41)./STATES(:,18));
    ALGEBRAIC(:,115) =  ( ( ( ( CONSTANTS(:,115).*CONSTANTS(:,53)).*power(STATES(:,34), 3.00000)).*STATES(:,32)).*STATES(:,33)).*(STATES(:,21) - ALGEBRAIC(:,114));
    ALGEBRAIC(:,117) =  ( CONSTANTS(:,115).*CONSTANTS(:,54)).*(STATES(:,21) - ALGEBRAIC(:,114));
    ALGEBRAIC(:,79) = (( ( ( CONSTANTS(:,115).*CONSTANTS(:,63)).*ALGEBRAIC(:,76)).*CONSTANTS(:,40))./(1.00000+power(CONSTANTS(:,65)./STATES(:,18), 4.00000)))./(CONSTANTS(:,40)+CONSTANTS(:,64));
    ALGEBRAIC(:,119) =  ( ( ( CONSTANTS(:,115).*CONSTANTS(:,107)).*power(STATES(:,36), 3.00000)).*STATES(:,35)).*(STATES(:,21) - ALGEBRAIC(:,114));
    ALGEBRAIC(:,123) = ((((ALGEBRAIC(:,115)+ALGEBRAIC(:,117))+ 3.00000.*ALGEBRAIC(:,75))+ 3.00000.*ALGEBRAIC(:,79))+ALGEBRAIC(:,49))+ALGEBRAIC(:,119);
    RATES(:,18) = (((  - ALGEBRAIC(:,123).*CONSTANTS(:,28))./( CONSTANTS(:,121).*CONSTANTS(:,29))+ (CONSTANTS(:,111)./CONSTANTS(:,121)).*(STATES(:,17) - STATES(:,18)))+ (CONSTANTS(:,112)./CONSTANTS(:,121)).*(STATES(:,28) - STATES(:,18))) - RATES(:,16);
    ALGEBRAIC(:,126) =  piecewise({(VOI - CONSTANTS(:,86)) -  CONSTANTS(:,87).*floor((VOI - CONSTANTS(:,86))./CONSTANTS(:,87))<CONSTANTS(:,85), 1.00000 }, 0.00000).*CONSTANTS(:,113);
    ALGEBRAIC(:,97) = ALGEBRAIC(:,89)+ALGEBRAIC(:,96);
    ALGEBRAIC(:,50) =  CONSTANTS(:,45).*(STATES(:,21) - CONSTANTS(:,116));
    ALGEBRAIC(:,52) =  (( CONSTANTS(:,36).*CONSTANTS(:,47))./(1.00000+CONSTANTS(:,46)./STATES(:,3))).*(STATES(:,21) - CONSTANTS(:,116));
    ALGEBRAIC(:,54) =  (( CONSTANTS(:,115).*CONSTANTS(:,47))./(1.00000+CONSTANTS(:,46)./STATES(:,4))).*(STATES(:,21) - CONSTANTS(:,116));
    ALGEBRAIC(:,55) = ALGEBRAIC(:,52)+ALGEBRAIC(:,54);
    ALGEBRAIC(:,85) = ALGEBRAIC(:,55)+ALGEBRAIC(:,50);
    ALGEBRAIC(:,44) = ( ( ( ( CONSTANTS(:,96).*STATES(:,21)).*CONSTANTS(:,29)).*CONSTANTS(:,93)).*( ( 0.750000.*STATES(:,22)).*exp( STATES(:,21).*CONSTANTS(:,93)) -  0.750000.*CONSTANTS(:,40)))./(exp( STATES(:,21).*CONSTANTS(:,93)) - 1.00000);
    ALGEBRAIC(:,45) =  ( ( ( ( ALGEBRAIC(:,44).*STATES(:,23)).*STATES(:,24)).*( CONSTANTS(:,39).*(CONSTANTS(:,44)+(1.00000 - STATES(:,25)))+ CONSTANTS(:,94).*(CONSTANTS(:,43)+(1.00000 - STATES(:,26))))).*power(CONSTANTS(:,42), CONSTANTS(:,118))).*0.450000;
    ALGEBRAIC(:,98) =  (1.00000./CONSTANTS(:,93)).*log(CONSTANTS(:,40)./STATES(:,22));
    ALGEBRAIC(:,56) = STATES(:,18);
    ALGEBRAIC(:,57) = 0.100000+0.900000./(1.00000+power(ALGEBRAIC(:,56)./7.00000, 2.00000));
    ALGEBRAIC(:,99) = ( ALGEBRAIC(:,57).*1.00000)./(1.00000+exp( 0.238500.*((STATES(:,21) - ALGEBRAIC(:,98)) - 59.2150)));
    ALGEBRAIC(:,100) =  1.00000.*(( 0.491240.*exp( 0.0803200.*((STATES(:,21) - ALGEBRAIC(:,98))+5.47600))+exp( 0.0617500.*((STATES(:,21) - ALGEBRAIC(:,98)) - 594.310)))./(1.00000+exp(  - 0.514300.*((STATES(:,21) - ALGEBRAIC(:,98))+4.75300))));
    ALGEBRAIC(:,101) = ALGEBRAIC(:,99)./(ALGEBRAIC(:,99)+ALGEBRAIC(:,100));
    ALGEBRAIC(:,102) =  ( ( CONSTANTS(:,48).*CONSTANTS(:,98)).*ALGEBRAIC(:,101)).*(STATES(:,21) - ALGEBRAIC(:,98));
    ALGEBRAIC(:,58) = STATES(:,18);
    ALGEBRAIC(:,59) = piecewise({CONSTANTS(:,32)<1.00000, 1.00000+1.50000./(1.00000+power(9.00000./ALGEBRAIC(:,58), 4.00000)) }, 1.00000);
    ALGEBRAIC(:,103) =  ( ( ( ( CONSTANTS(:,100).*ALGEBRAIC(:,59)).*CONSTANTS(:,99)).*(0.0550000+0.400000./(1.00000+exp(((STATES(:,21) - ALGEBRAIC(:,98))+9.53000)./17.1800)))).*(STATES(:,21) - ALGEBRAIC(:,98))).*power((CONSTANTS(:,40)./5.40000), 1.0 ./ 2);
    ALGEBRAIC(:,60) = 1.00000./(1.00000+exp(7.48800 - STATES(:,21)./5.98000));
    ALGEBRAIC(:,104) =  ( ( CONSTANTS(:,36).*CONSTANTS(:,51)).*ALGEBRAIC(:,60)).*(STATES(:,21) - ALGEBRAIC(:,98));
    ALGEBRAIC(:,105) =  ( ( CONSTANTS(:,115).*CONSTANTS(:,51)).*ALGEBRAIC(:,60)).*(STATES(:,21) - ALGEBRAIC(:,98));
    ALGEBRAIC(:,106) = ALGEBRAIC(:,104)+ALGEBRAIC(:,105);
    ALGEBRAIC(:,61) = 1.00000./(1.00000+exp((STATES(:,21)+74.0000)./24.0000));
    ALGEBRAIC(:,107) =  ( ( CONSTANTS(:,101).*STATES(:,27)).*ALGEBRAIC(:,61)).*(STATES(:,21) - ALGEBRAIC(:,98));
    ALGEBRAIC(:,62) =  (1.00000./CONSTANTS(:,93)).*log((CONSTANTS(:,40)+ CONSTANTS(:,52).*CONSTANTS(:,41))./(STATES(:,22)+ CONSTANTS(:,52).*STATES(:,28)));
    ALGEBRAIC(:,63) =  ( ( CONSTANTS(:,36).*CONSTANTS(:,102)).*power(STATES(:,29), 2.00000)).*(STATES(:,21) - ALGEBRAIC(:,62));
    ALGEBRAIC(:,64) =  ( ( CONSTANTS(:,115).*CONSTANTS(:,103)).*power(STATES(:,29), 2.00000)).*(STATES(:,21) - ALGEBRAIC(:,62));
    ALGEBRAIC(:,65) = ALGEBRAIC(:,63)+ALGEBRAIC(:,64);
    ALGEBRAIC(:,108) =  ( ( CONSTANTS(:,104).*STATES(:,30)).*STATES(:,31)).*(STATES(:,21) - ALGEBRAIC(:,98));
    ALGEBRAIC(:,80) = ALGEBRAIC(:,78)+ALGEBRAIC(:,79);
    ALGEBRAIC(:,109) =  ( ( CONSTANTS(:,108).*STATES(:,37)).*STATES(:,38)).*(STATES(:,21) - ALGEBRAIC(:,98));
    ALGEBRAIC(:,121) = (((((((ALGEBRAIC(:,109)+ALGEBRAIC(:,107))+ALGEBRAIC(:,65))+ALGEBRAIC(:,102)) -  2.00000.*ALGEBRAIC(:,80))+ALGEBRAIC(:,45))+ALGEBRAIC(:,106))+ALGEBRAIC(:,108))+ALGEBRAIC(:,103);
    ALGEBRAIC(:,124) = ALGEBRAIC(:,122)+ALGEBRAIC(:,123);
    ALGEBRAIC(:,125) = ((ALGEBRAIC(:,124)+ALGEBRAIC(:,85))+ALGEBRAIC(:,97))+ALGEBRAIC(:,121);
    RATES(:,21) =  - (ALGEBRAIC(:,125)+ALGEBRAIC(:,126));
   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(:,10) = 1.00000./(1.00000+exp((STATES(:,21)+91.0000)./6.10000));
    ALGEBRAIC(:,1) = 1.00000./(1.00000+exp( - (STATES(:,21)+9.00000)./6.00000));
    ALGEBRAIC(:,15) = ( ALGEBRAIC(:,1).*(1.00000 - exp( - (STATES(:,21)+9.00000)./6.00000)))./( 0.0350000.*(STATES(:,21)+9.00000));
    ALGEBRAIC(:,2) = 1.00000./(1.00000+exp((STATES(:,21)+30.0000)./7.00000))+0.200000./(1.00000+exp((50.0000 - STATES(:,21))./20.0000));
    ALGEBRAIC(:,16) = 1.00000./( 0.0197000.*exp( - power( 0.0337000.*(STATES(:,21)+25.0000), 2.00000))+0.0200000);
    ALGEBRAIC(:,3) = 1.00000./(1.00000+exp( - (STATES(:,21)+10.0000)./5.00000));
    ALGEBRAIC(:,17) = ( ((550.000)./(1.00000+exp(( - 22.0000 - STATES(:,21))./9.00000))).*6.00000)./(1.00000+exp((STATES(:,21)+11.0000)./9.00000))+230.000./(1.00000+exp((STATES(:,21)+40.0000)./20.0000));
    ALGEBRAIC(:,4) = 1.00000./(1.00000+exp( - (STATES(:,21)+3.80000)./14.2500));
    ALGEBRAIC(:,18) = 990.100./(1.00000+exp( - (STATES(:,21)+2.43600)./14.1200));
    ALGEBRAIC(:,5) = 1.00000./(1.00000+exp((STATES(:,21)+6.00000)./ - 8.60000));
    ALGEBRAIC(:,19) = 9.00000./(1.00000+exp((STATES(:,21)+5.00000)./12.0000))+0.500000;
    ALGEBRAIC(:,6) = 1.00000./(1.00000+exp((STATES(:,21)+7.50000)./10.0000));
    ALGEBRAIC(:,20) = 590.000./(1.00000+exp((STATES(:,21)+60.0000)./10.0000))+3050.00;
    ALGEBRAIC(:,11) = ( 0.320000.*(STATES(:,21)+47.1300))./(1.00000 - exp(  - 0.100000.*(STATES(:,21)+47.1300)));
    ALGEBRAIC(:,24) =  0.0800000.*exp( - STATES(:,21)./11.0000);
    ALGEBRAIC(:,12) = 1.00000./(1.00000+exp( - (STATES(:,21)+1.00000)./11.0000));
    ALGEBRAIC(:,25) =  3.50000.*exp( - power(STATES(:,21)./30.0000, 2.00000))+1.50000;
    ALGEBRAIC(:,13) = 1.00000./(1.00000+exp((STATES(:,21)+40.5000)./11.5000));
    ALGEBRAIC(:,26) =  25.6350.*exp( - power((STATES(:,21)+52.4500)./15.8827, 2.00000))+24.1400;
    ALGEBRAIC(:,7) = piecewise({STATES(:,21)< - 40.0000,  0.135000.*exp((STATES(:,21)+80.0000)./ - 6.80000) }, 0.00000);
    ALGEBRAIC(:,21) = piecewise({STATES(:,21)< - 40.0000,  3.56000.*exp( 0.0790000.*STATES(:,21))+ 310000..*exp( 0.350000.*STATES(:,21)) }, 1.00000./( 0.130000.*(1.00000+exp((STATES(:,21)+10.6600)./ - 11.1000))));
    ALGEBRAIC(:,28) = ALGEBRAIC(:,7)./(ALGEBRAIC(:,7)+ALGEBRAIC(:,21));
    ALGEBRAIC(:,32) = 1.00000./(ALGEBRAIC(:,7)+ALGEBRAIC(:,21));
    ALGEBRAIC(:,8) = piecewise({STATES(:,21)< - 40.0000, ( (  - 127140..*exp( 0.244400.*STATES(:,21)) -  3.47400e-05.*exp(  - 0.0439100.*STATES(:,21))).*(STATES(:,21)+37.7800))./(1.00000+exp( 0.311000.*(STATES(:,21)+79.2300))) }, 0.00000);
    ALGEBRAIC(:,22) = piecewise({STATES(:,21)< - 40.0000, ( 0.121200.*exp(  - 0.0105200.*STATES(:,21)))./(1.00000+exp(  - 0.137800.*(STATES(:,21)+40.1400))) }, ( 0.300000.*exp(  - 2.53500e-07.*STATES(:,21)))./(1.00000+exp(  - 0.100000.*(STATES(:,21)+32.0000))));
    ALGEBRAIC(:,29) = ALGEBRAIC(:,8)./(ALGEBRAIC(:,8)+ALGEBRAIC(:,22));
    ALGEBRAIC(:,33) = 1.00000./(ALGEBRAIC(:,8)+ALGEBRAIC(:,22));
    ALGEBRAIC(:,9) = piecewise({STATES(:,21)== - 47.1300, 3.20000 }, ( 0.320000.*(STATES(:,21)+47.1300))./(1.00000 - exp(  - 0.100000.*(STATES(:,21)+47.1300))));
    ALGEBRAIC(:,23) =  0.0800000.*exp( - STATES(:,21)./11.0000);
    ALGEBRAIC(:,30) = ALGEBRAIC(:,9)./(ALGEBRAIC(:,9)+ALGEBRAIC(:,23));
    ALGEBRAIC(:,34) = 1.00000./(ALGEBRAIC(:,9)+ALGEBRAIC(:,23));
    ALGEBRAIC(:,14) = ((1.00000 - STATES(:,41)) - STATES(:,40)) - STATES(:,39);
    ALGEBRAIC(:,27) = CONSTANTS(:,75) - (CONSTANTS(:,75) - CONSTANTS(:,76))./(1.00000+power(CONSTANTS(:,79)./STATES(:,19), 2.50000));
    ALGEBRAIC(:,31) =  CONSTANTS(:,81).*ALGEBRAIC(:,27);
    ALGEBRAIC(:,35) = CONSTANTS(:,110)./ALGEBRAIC(:,27);
    ALGEBRAIC(:,36) = RATES(:,9)+RATES(:,7);
    ALGEBRAIC(:,87) =  ((1.00000./CONSTANTS(:,93))./2.00000).*log(CONSTANTS(:,38)./STATES(:,3));
    ALGEBRAIC(:,88) =  ( CONSTANTS(:,36).*CONSTANTS(:,37)).*(STATES(:,21) - ALGEBRAIC(:,87));
    ALGEBRAIC(:,39) = ( ( ( ( ( CONSTANTS(:,95).*4.00000).*STATES(:,21)).*CONSTANTS(:,29)).*CONSTANTS(:,93)).*( ( 0.341000.*STATES(:,3)).*exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,93)) -  0.341000.*CONSTANTS(:,38)))./(exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,93)) - 1.00000);
    ALGEBRAIC(:,40) =  ( ( ( ( ( CONSTANTS(:,39).*ALGEBRAIC(:,39)).*STATES(:,23)).*STATES(:,24)).*((1.00000 - STATES(:,25))+CONSTANTS(:,44))).*power(CONSTANTS(:,42), CONSTANTS(:,118))).*0.450000;
    ALGEBRAIC(:,66) = 1.00000./(1.00000+power(CONSTANTS(:,55)./STATES(:,3), 2.00000));
    ALGEBRAIC(:,68) =  ( exp( ( CONSTANTS(:,62).*STATES(:,21)).*CONSTANTS(:,93)).*power(STATES(:,17), 3.00000)).*CONSTANTS(:,38);
    ALGEBRAIC(:,69) =  ( exp( ( (CONSTANTS(:,62) - 1.00000).*STATES(:,21)).*CONSTANTS(:,93)).*power(CONSTANTS(:,41), 3.00000)).*STATES(:,3);
    ALGEBRAIC(:,70) = ((( ( CONSTANTS(:,56).*power(CONSTANTS(:,41), 3.00000)).*(1.00000+power(STATES(:,17)./CONSTANTS(:,58), 3.00000))+ ( power(CONSTANTS(:,59), 3.00000).*STATES(:,3)).*(1.00000+STATES(:,3)./CONSTANTS(:,56)))+ CONSTANTS(:,57).*power(STATES(:,17), 3.00000))+ power(STATES(:,17), 3.00000).*CONSTANTS(:,38))+ power(CONSTANTS(:,41), 3.00000).*STATES(:,3);
    ALGEBRAIC(:,71) = (( ( ( ( CONSTANTS(:,36).*CONSTANTS(:,105)).*power(CONSTANTS(:,60), CONSTANTS(:,118))).*ALGEBRAIC(:,66)).*(ALGEBRAIC(:,68) - ALGEBRAIC(:,69)))./ALGEBRAIC(:,70))./(1.00000+ CONSTANTS(:,61).*exp( ( (CONSTANTS(:,62) - 1.00000).*STATES(:,21)).*CONSTANTS(:,93)));
    ALGEBRAIC(:,81) = power( STATES(:,3).*1.00000, 1.60000);
    ALGEBRAIC(:,82) = ( ( ( CONSTANTS(:,36).*power(CONSTANTS(:,71), CONSTANTS(:,118))).*CONSTANTS(:,69)).*ALGEBRAIC(:,81))./(ALGEBRAIC(:,81)+CONSTANTS(:,90));
    ALGEBRAIC(:,89) = ((ALGEBRAIC(:,40)+ALGEBRAIC(:,88))+ALGEBRAIC(:,82)) -  2.00000.*ALGEBRAIC(:,71);
    ALGEBRAIC(:,90) =  ( CONSTANTS(:,84).*STATES(:,40)).*(STATES(:,19) - STATES(:,3));
    ALGEBRAIC(:,92) =  ( (1.00000+ 0.250000.*CONSTANTS(:,32)).*(STATES(:,19) - STATES(:,3))).*5.34800e-06;
    ALGEBRAIC(:,38) = (((((RATES(:,14)+RATES(:,12))+RATES(:,13))+RATES(:,1))+RATES(:,5))+RATES(:,6))+RATES(:,11);
    ALGEBRAIC(:,94) = ( ( power(CONSTANTS(:,77), CONSTANTS(:,118)).*CONSTANTS(:,78)).*(power(STATES(:,2)./CONSTANTS(:,109), CONSTANTS(:,80)) - power(STATES(:,19)./CONSTANTS(:,74), CONSTANTS(:,80))))./((1.00000+power(STATES(:,2)./CONSTANTS(:,109), CONSTANTS(:,80)))+power(STATES(:,19)./CONSTANTS(:,74), CONSTANTS(:,80)));
    ALGEBRAIC(:,37) = RATES(:,10)+RATES(:,8);
    ALGEBRAIC(:,91) =  ((1.00000./CONSTANTS(:,93))./2.00000).*log(CONSTANTS(:,38)./STATES(:,4));
    ALGEBRAIC(:,93) =  ( CONSTANTS(:,115).*CONSTANTS(:,37)).*(STATES(:,21) - ALGEBRAIC(:,91));
    ALGEBRAIC(:,41) = ( ( ( ( ( CONSTANTS(:,95).*4.00000).*STATES(:,21)).*CONSTANTS(:,29)).*CONSTANTS(:,93)).*( ( 0.341000.*STATES(:,4)).*exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,93)) -  0.341000.*CONSTANTS(:,38)))./(exp( ( 2.00000.*STATES(:,21)).*CONSTANTS(:,93)) - 1.00000);
    ALGEBRAIC(:,42) =  ( ( ( ( ( CONSTANTS(:,94).*ALGEBRAIC(:,41)).*STATES(:,23)).*STATES(:,24)).*((1.00000 - STATES(:,26))+CONSTANTS(:,43))).*power(CONSTANTS(:,42), CONSTANTS(:,118))).*0.450000;
    ALGEBRAIC(:,67) = 1.00000./(1.00000+power(CONSTANTS(:,55)./STATES(:,4), 2.00000));
    ALGEBRAIC(:,72) =  ( exp( ( CONSTANTS(:,62).*STATES(:,21)).*CONSTANTS(:,93)).*power(STATES(:,18), 3.00000)).*CONSTANTS(:,38);
    ALGEBRAIC(:,73) =  ( exp( ( (CONSTANTS(:,62) - 1.00000).*STATES(:,21)).*CONSTANTS(:,93)).*power(CONSTANTS(:,41), 3.00000)).*STATES(:,4);
    ALGEBRAIC(:,74) = ((( ( CONSTANTS(:,56).*power(CONSTANTS(:,41), 3.00000)).*(1.00000+power(STATES(:,18)./CONSTANTS(:,58), 3.00000))+ ( power(CONSTANTS(:,59), 3.00000).*STATES(:,4)).*(1.00000+STATES(:,4)./CONSTANTS(:,56)))+ CONSTANTS(:,57).*power(STATES(:,18), 3.00000))+ power(STATES(:,18), 3.00000).*CONSTANTS(:,38))+ power(CONSTANTS(:,41), 3.00000).*STATES(:,4);
    ALGEBRAIC(:,75) = (( ( ( ( CONSTANTS(:,115).*CONSTANTS(:,105)).*power(CONSTANTS(:,60), CONSTANTS(:,118))).*ALGEBRAIC(:,67)).*(ALGEBRAIC(:,72) - ALGEBRAIC(:,73)))./ALGEBRAIC(:,74))./(1.00000+ CONSTANTS(:,61).*exp( ( (CONSTANTS(:,62) - 1.00000).*STATES(:,21)).*CONSTANTS(:,93)));
    ALGEBRAIC(:,83) = power( STATES(:,4).*1.00000, 1.60000);
    ALGEBRAIC(:,84) = ( ( ( CONSTANTS(:,115).*power(CONSTANTS(:,71), CONSTANTS(:,118))).*CONSTANTS(:,69)).*ALGEBRAIC(:,83))./(ALGEBRAIC(:,83)+CONSTANTS(:,90));
    ALGEBRAIC(:,96) = ((ALGEBRAIC(:,42)+ALGEBRAIC(:,93))+ALGEBRAIC(:,84)) -  2.00000.*ALGEBRAIC(:,75);
    ALGEBRAIC(:,46) = ( ( ( ( CONSTANTS(:,97).*STATES(:,21)).*CONSTANTS(:,29)).*CONSTANTS(:,93)).*( ( 0.750000.*STATES(:,17)).*exp( STATES(:,21).*CONSTANTS(:,93)) -  0.750000.*CONSTANTS(:,41)))./(exp( STATES(:,21).*CONSTANTS(:,93)) - 1.00000);
    ALGEBRAIC(:,47) =  ( ( ( ( ( CONSTANTS(:,39).*ALGEBRAIC(:,46)).*STATES(:,23)).*STATES(:,24)).*((1.00000 - STATES(:,25))+CONSTANTS(:,44))).*power(CONSTANTS(:,42), CONSTANTS(:,118))).*0.450000;
    ALGEBRAIC(:,110) =  (1.00000./CONSTANTS(:,93)).*log(CONSTANTS(:,41)./STATES(:,17));
    ALGEBRAIC(:,111) =  ( ( ( ( CONSTANTS(:,36).*CONSTANTS(:,53)).*power(STATES(:,34), 3.00000)).*STATES(:,32)).*STATES(:,33)).*(STATES(:,21) - ALGEBRAIC(:,110));
    ALGEBRAIC(:,112) =  ( CONSTANTS(:,36).*CONSTANTS(:,54)).*(STATES(:,21) - ALGEBRAIC(:,110));
    ALGEBRAIC(:,76) = 1.00000./((1.00000+ 0.124500.*exp( (  - 0.100000.*STATES(:,21)).*CONSTANTS(:,93)))+ ( 0.0365000.*CONSTANTS(:,106)).*exp(  - STATES(:,21).*CONSTANTS(:,93)));
    ALGEBRAIC(:,78) = (( ( ( CONSTANTS(:,36).*CONSTANTS(:,63)).*ALGEBRAIC(:,76)).*CONSTANTS(:,40))./(1.00000+power(CONSTANTS(:,65)./STATES(:,17), 4.00000)))./(CONSTANTS(:,40)+CONSTANTS(:,64));
    ALGEBRAIC(:,113) =  ( ( ( CONSTANTS(:,36).*CONSTANTS(:,107)).*power(STATES(:,36), 3.00000)).*STATES(:,35)).*(STATES(:,21) - ALGEBRAIC(:,110));
    ALGEBRAIC(:,122) = ((((ALGEBRAIC(:,111)+ALGEBRAIC(:,112))+ 3.00000.*ALGEBRAIC(:,71))+ 3.00000.*ALGEBRAIC(:,78))+ALGEBRAIC(:,47))+ALGEBRAIC(:,113);
    ALGEBRAIC(:,48) = ( ( ( ( CONSTANTS(:,97).*STATES(:,21)).*CONSTANTS(:,29)).*CONSTANTS(:,93)).*( ( 0.750000.*STATES(:,18)).*exp( STATES(:,21).*CONSTANTS(:,93)) -  0.750000.*CONSTANTS(:,41)))./(exp( STATES(:,21).*CONSTANTS(:,93)) - 1.00000);
    ALGEBRAIC(:,49) =  ( ( ( ( ( CONSTANTS(:,94).*ALGEBRAIC(:,48)).*STATES(:,23)).*STATES(:,24)).*((1.00000 - STATES(:,26))+CONSTANTS(:,43))).*power(CONSTANTS(:,42), CONSTANTS(:,118))).*0.450000;
    ALGEBRAIC(:,114) =  (1.00000./CONSTANTS(:,93)).*log(CONSTANTS(:,41)./STATES(:,18));
    ALGEBRAIC(:,115) =  ( ( ( ( CONSTANTS(:,115).*CONSTANTS(:,53)).*power(STATES(:,34), 3.00000)).*STATES(:,32)).*STATES(:,33)).*(STATES(:,21) - ALGEBRAIC(:,114));
    ALGEBRAIC(:,117) =  ( CONSTANTS(:,115).*CONSTANTS(:,54)).*(STATES(:,21) - ALGEBRAIC(:,114));
    ALGEBRAIC(:,79) = (( ( ( CONSTANTS(:,115).*CONSTANTS(:,63)).*ALGEBRAIC(:,76)).*CONSTANTS(:,40))./(1.00000+power(CONSTANTS(:,65)./STATES(:,18), 4.00000)))./(CONSTANTS(:,40)+CONSTANTS(:,64));
    ALGEBRAIC(:,119) =  ( ( ( CONSTANTS(:,115).*CONSTANTS(:,107)).*power(STATES(:,36), 3.00000)).*STATES(:,35)).*(STATES(:,21) - ALGEBRAIC(:,114));
    ALGEBRAIC(:,123) = ((((ALGEBRAIC(:,115)+ALGEBRAIC(:,117))+ 3.00000.*ALGEBRAIC(:,75))+ 3.00000.*ALGEBRAIC(:,79))+ALGEBRAIC(:,49))+ALGEBRAIC(:,119);
    ALGEBRAIC(:,126) =  piecewise({(VOI - CONSTANTS(:,86)) -  CONSTANTS(:,87).*floor((VOI - CONSTANTS(:,86))./CONSTANTS(:,87))<CONSTANTS(:,85), 1.00000 }, 0.00000).*CONSTANTS(:,113);
    ALGEBRAIC(:,97) = ALGEBRAIC(:,89)+ALGEBRAIC(:,96);
    ALGEBRAIC(:,50) =  CONSTANTS(:,45).*(STATES(:,21) - CONSTANTS(:,116));
    ALGEBRAIC(:,52) =  (( CONSTANTS(:,36).*CONSTANTS(:,47))./(1.00000+CONSTANTS(:,46)./STATES(:,3))).*(STATES(:,21) - CONSTANTS(:,116));
    ALGEBRAIC(:,54) =  (( CONSTANTS(:,115).*CONSTANTS(:,47))./(1.00000+CONSTANTS(:,46)./STATES(:,4))).*(STATES(:,21) - CONSTANTS(:,116));
    ALGEBRAIC(:,55) = ALGEBRAIC(:,52)+ALGEBRAIC(:,54);
    ALGEBRAIC(:,85) = ALGEBRAIC(:,55)+ALGEBRAIC(:,50);
    ALGEBRAIC(:,44) = ( ( ( ( CONSTANTS(:,96).*STATES(:,21)).*CONSTANTS(:,29)).*CONSTANTS(:,93)).*( ( 0.750000.*STATES(:,22)).*exp( STATES(:,21).*CONSTANTS(:,93)) -  0.750000.*CONSTANTS(:,40)))./(exp( STATES(:,21).*CONSTANTS(:,93)) - 1.00000);
    ALGEBRAIC(:,45) =  ( ( ( ( ALGEBRAIC(:,44).*STATES(:,23)).*STATES(:,24)).*( CONSTANTS(:,39).*(CONSTANTS(:,44)+(1.00000 - STATES(:,25)))+ CONSTANTS(:,94).*(CONSTANTS(:,43)+(1.00000 - STATES(:,26))))).*power(CONSTANTS(:,42), CONSTANTS(:,118))).*0.450000;
    ALGEBRAIC(:,98) =  (1.00000./CONSTANTS(:,93)).*log(CONSTANTS(:,40)./STATES(:,22));
    ALGEBRAIC(:,56) = STATES(:,18);
    ALGEBRAIC(:,57) = 0.100000+0.900000./(1.00000+power(ALGEBRAIC(:,56)./7.00000, 2.00000));
    ALGEBRAIC(:,99) = ( ALGEBRAIC(:,57).*1.00000)./(1.00000+exp( 0.238500.*((STATES(:,21) - ALGEBRAIC(:,98)) - 59.2150)));
    ALGEBRAIC(:,100) =  1.00000.*(( 0.491240.*exp( 0.0803200.*((STATES(:,21) - ALGEBRAIC(:,98))+5.47600))+exp( 0.0617500.*((STATES(:,21) - ALGEBRAIC(:,98)) - 594.310)))./(1.00000+exp(  - 0.514300.*((STATES(:,21) - ALGEBRAIC(:,98))+4.75300))));
    ALGEBRAIC(:,101) = ALGEBRAIC(:,99)./(ALGEBRAIC(:,99)+ALGEBRAIC(:,100));
    ALGEBRAIC(:,102) =  ( ( CONSTANTS(:,48).*CONSTANTS(:,98)).*ALGEBRAIC(:,101)).*(STATES(:,21) - ALGEBRAIC(:,98));
    ALGEBRAIC(:,58) = STATES(:,18);
    ALGEBRAIC(:,59) = piecewise({CONSTANTS(:,32)<1.00000, 1.00000+1.50000./(1.00000+power(9.00000./ALGEBRAIC(:,58), 4.00000)) }, 1.00000);
    ALGEBRAIC(:,103) =  ( ( ( ( CONSTANTS(:,100).*ALGEBRAIC(:,59)).*CONSTANTS(:,99)).*(0.0550000+0.400000./(1.00000+exp(((STATES(:,21) - ALGEBRAIC(:,98))+9.53000)./17.1800)))).*(STATES(:,21) - ALGEBRAIC(:,98))).*power((CONSTANTS(:,40)./5.40000), 1.0 ./ 2);
    ALGEBRAIC(:,60) = 1.00000./(1.00000+exp(7.48800 - STATES(:,21)./5.98000));
    ALGEBRAIC(:,104) =  ( ( CONSTANTS(:,36).*CONSTANTS(:,51)).*ALGEBRAIC(:,60)).*(STATES(:,21) - ALGEBRAIC(:,98));
    ALGEBRAIC(:,105) =  ( ( CONSTANTS(:,115).*CONSTANTS(:,51)).*ALGEBRAIC(:,60)).*(STATES(:,21) - ALGEBRAIC(:,98));
    ALGEBRAIC(:,106) = ALGEBRAIC(:,104)+ALGEBRAIC(:,105);
    ALGEBRAIC(:,61) = 1.00000./(1.00000+exp((STATES(:,21)+74.0000)./24.0000));
    ALGEBRAIC(:,107) =  ( ( CONSTANTS(:,101).*STATES(:,27)).*ALGEBRAIC(:,61)).*(STATES(:,21) - ALGEBRAIC(:,98));
    ALGEBRAIC(:,62) =  (1.00000./CONSTANTS(:,93)).*log((CONSTANTS(:,40)+ CONSTANTS(:,52).*CONSTANTS(:,41))./(STATES(:,22)+ CONSTANTS(:,52).*STATES(:,28)));
    ALGEBRAIC(:,63) =  ( ( CONSTANTS(:,36).*CONSTANTS(:,102)).*power(STATES(:,29), 2.00000)).*(STATES(:,21) - ALGEBRAIC(:,62));
    ALGEBRAIC(:,64) =  ( ( CONSTANTS(:,115).*CONSTANTS(:,103)).*power(STATES(:,29), 2.00000)).*(STATES(:,21) - ALGEBRAIC(:,62));
    ALGEBRAIC(:,65) = ALGEBRAIC(:,63)+ALGEBRAIC(:,64);
    ALGEBRAIC(:,108) =  ( ( CONSTANTS(:,104).*STATES(:,30)).*STATES(:,31)).*(STATES(:,21) - ALGEBRAIC(:,98));
    ALGEBRAIC(:,80) = ALGEBRAIC(:,78)+ALGEBRAIC(:,79);
    ALGEBRAIC(:,109) =  ( ( CONSTANTS(:,108).*STATES(:,37)).*STATES(:,38)).*(STATES(:,21) - ALGEBRAIC(:,98));
    ALGEBRAIC(:,121) = (((((((ALGEBRAIC(:,109)+ALGEBRAIC(:,107))+ALGEBRAIC(:,65))+ALGEBRAIC(:,102)) -  2.00000.*ALGEBRAIC(:,80))+ALGEBRAIC(:,45))+ALGEBRAIC(:,106))+ALGEBRAIC(:,108))+ALGEBRAIC(:,103);
    ALGEBRAIC(:,124) = ALGEBRAIC(:,122)+ALGEBRAIC(:,123);
    ALGEBRAIC(:,125) = ((ALGEBRAIC(:,124)+ALGEBRAIC(:,85))+ALGEBRAIC(:,97))+ALGEBRAIC(:,121);
    ALGEBRAIC(:,43) = ALGEBRAIC(:,40)+ALGEBRAIC(:,42);
    ALGEBRAIC(:,51) = ALGEBRAIC(:,47)+ALGEBRAIC(:,49);
    ALGEBRAIC(:,53) = (ALGEBRAIC(:,43)+ALGEBRAIC(:,51))+ALGEBRAIC(:,45);
    ALGEBRAIC(:,77) = ALGEBRAIC(:,71)+ALGEBRAIC(:,75);
    ALGEBRAIC(:,86) = ALGEBRAIC(:,82)+ALGEBRAIC(:,84);
    ALGEBRAIC(:,95) = ALGEBRAIC(:,88)+ALGEBRAIC(:,93);
    ALGEBRAIC(:,116) = ALGEBRAIC(:,111)+ALGEBRAIC(:,115);
    ALGEBRAIC(:,118) = ALGEBRAIC(:,112)+ALGEBRAIC(:,117);
    ALGEBRAIC(:,120) = ALGEBRAIC(:,113)+ALGEBRAIC(:,119);
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

% 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